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ABSTRACT 


The' Multiple Model Adaptive Control algorithm has been used m appli- 
cations of advanced control technology. However, in these applications, 
many undesirable characteristics of the method, such as high amplitude 
limit cycles, have been uncovered. In this thesis the basic types 
of behavior exhibited by the method are explored. This is done through 
the simulation and analysis of the method as applied to a sample system 
structure. This structure has been carefully chosen to exhibit the 
major phenomena of interest while remaining amenable to detailed analysis. 
Two major types of results are presented. First of all, detailed conditions 
for the existance of each of the types of behavior are developed for the 
special system structure under consideration. Of possibly greater signi- 
ficance are the qualitative insights which result from extrapolating the 
detailed conclusions to problems of more general structure. It is believed 
that the qualitative understandings developed in this thesis can form the 
basis for the introduction of design modifications (two of which are 
suggested m this thesis) and the development of a systematic methodology 
for the design of adaptive control systems using the Multiple Model 
Adaptive Control algorithms. 


THESIS SUPERVISOR: Alan S. Willsky 

TITLE: Associate Professor, Electrical Engineering and 
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CHAPTER 1 


INTRODUCTION 


1. 1 Motivation 

In many applications of control theory, the dynamics of the plant 
are incompletely known at best. Furthermore, the dynamics are often time- 
varying and non-linear. In such an environment, control becomes a very 
difficult task and the problem of the optimal control of such systems 
remains unsolved. However, such systems need to be controlled and so 
a myriad of suboptimal schemes have emerged. 

The many methods which have been proposed can, in general, be di- 
vided into too classes; the passive methods which rely on the robustness 
of a time- invariant feedback controller to maintain good performance and 
the active methods which involve changing the controllers as necessary. 

For example, much of the work of Horowitz [1,2] has been aimed at de- 
riving a single, time-invariant control law which gives acceptable be- 
havior for all plant parameter values (a passive approach) . The work 
of Wong [3, 4] has similarly been aimed at analyzing the robustness 
properties of feedback controllers using a geometric approach, and 
Safonov [5] has derived robustness conditions for controllers when the 
parameter variations are due to a change in the operating point of a non- 
linear system. It should be pointed out that the ad hoc approach of 
increasing the plant noise design parameter (see Section 2.1) often 
mentioned for the standard Linear-Quadratic-Gauss ian (LQG) problem [6] 
is also a passive method of overcoming plant uncertainty. 

A major problem with such methods is that they are, by design, 
compromises. Performance for normal conditions is sacrificed in order 
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to improve performance for other conditions. In the extreme/ it may be 
impossible to maintain the desired performance for the full range of con- 
ditions" using - a fixed controller. 

In contrast to the passive methods, the active methods make use of 
a time- varying controller. Thus, they employ mechanisms which force the 
controller to adapt to changes m the operating environment. This, at 
least in theory, improves performance under all operating conditions 
since the controller can be tuned to the actual , rather than the average 
or even worst case, conditions. 

As previously mentioned, the optimal control of such systems is an 
unsolved problem. Thus ad hoc, suboptimal techniques have been proposed. 
However, partly because of the non-linearities (due to adaption) of such 
methods, they have not been subject to careful study regarding qualitative 
performance characteristics such as deterministic stability. In applica- 
tions, many of these methods have exhibited difficulties which have been 
mitigated by further ad hoc modifications of the design [16, 232 - In 
general, these modifications were not the product of extensive, systematic 
analysis of the system' s behavior and no_ general design methodology has 
emerged. 

The research which, is reported herein attempts to qualitatively and 
quantitatively examine the properties of one method of adaptive control 
which has been discussed in the literature, namely, the Multiple Model 
Adaptive Control (MMAC) method £7J . The MMAC method, which is discussed 
further in Chapter 2 of this thesis, has a very pleasing structure: a 
cascade of something which resembles a Maximum Aposteriori Probability 
(MAPI identifier I15J (basically a bank of Kalman Filters) and a bank of 
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linear quadratic regulators. However, in use it has become clear that 
the MMAC method can exhibit unacceptable behavior - such as high amplitude 
limit cycles. In this thesis the major qualitative properties of the 
MMAC method are examined and the principle reasons for the unacceptable 
behavior explored. It is believed that a through understanding of the 
behavior will lead to guidelines for the modification of the design 
which will ultimately yield a general design methodology. 

1.2 Background 

The general area of adaptive control has received much attention 
recently. For example [26J and [27.] both contain numerous references 
to a wide variety of approaches. The basic subject which is addressed 
is to generate a control uttl for a system given by 

xCt) = A(_t)xlt) + B(t)u{t) + £(.t) 
with observations 

y_(t] = C_tt)x(t) + nCt) 

The state xCtl is an n-vector while the input u(t) is an m-vector and 
y Ct) is a p-vector. The vectors (x) and rj_(t) represent system un- 
certainties and observation noise respectively. The condition which 
introduces the most complexity is that the system matrix A (nxn) , input 
matrix B (nxm) and output matrix C (pxn) are only incompletely known. 

The performance measure which is often used no judge such systems is 


a quadratic one, given by 
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ra 

J(u) = / [x‘ (t)Q_ x(t) + u' (t)Ru (t)]dt 

0 

where £ is an nxn -positive semi-defini-te matrix and- R is a mxm positive 
definite matrix. 

The solution to this control problem has not been found and the 
work in [7] clearly indicates that the optimal system, assuming it can 
be found, will prove to be far too complex to implement. Thus, numerous 
suboptimal solutions have been proposed. 

The MMAC method was largely inspired by the work of Magill [8] who 
was the first to examine a parallel -adaptive estimation algorithm in 
which the basic estimation is done by a bank of filters which are then 
coordinated by a centralized controller (see Figure 1.1). Further work 
in the area has been done by Lainiotis whose work has been summarized 
in [9] . The major thrust of this work has been aimed at adaptive esti- 
mation and parameter identification and not the control problem. 

Many authors [10, 11, 7] have examined feedback controllers based 
on the structure of Magill' s estimator, however. For example, Stein 
[10] has been able to derive upper and lower bounds on the cost m the 
optimal control problem and using the upper bound, has obtained a control 
law exhibiting a parallel structure. Saridis and Dao [11] have ex- 
ploited Stein’s lower bound to obtain a different control law. One 
major drawback of both algorithms is that they require significant 
on-line computation. 

Willner [7] proposed the MMAC algorithm as discussed m this 
thesis and showed that it performed well m relation to both the 




Fig. 1.1 Parallel-; 
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upper and lower bounds of Stein. Independently, Deshpande et al [12] 
arrived at the same algorithm as an ad hoc extension of the estimation/ 
identification algorithm of Mag'ill [8] . The method has been further 
considered by Lamiotis [13] . 

To our knowledge, no one has been successful in establishing any 
definitive properties on the behavior of the MMAC method. Specifically, 
the convergence properties of the identification problem (i.e., with 
the adaptive mechanism disabled) have only recently been proven- Both 
Hawkes and Moore [14] and Baram [15] have provided useful results, but 
both results do not hold in an adaptive situation (i.e., when the con- 
trol law is a function of the system output) . 

The MMAC method has been applied to various settings. For example, 
m Athans et al. [16, 23] the method has been used to control the F-8 
aircraft. Also, the F-8 controller of Stem et al. [17] can also be 
considered to be a multiple model design. Additional experience with 
and insight into the MMAC method has been gained by applying the esti- 
mation/identification algorithm to the detection of accidents on free- 

\ 

ways [18] . 

To a great extent, the F-8 project [16, 23] has provided the moti- 
vation for our work. In that project, where the true system does not 
correspond exactly to any of the hypothesized models, several problems 
were encountered. First of all, the probabilities often oscillated 
rapidly between two models, exhibiting behavior very much like a limit 
cycle. This problem led to a design modification consisting of the 
insertion of low pass filters to slow the probability transitions. 
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A second problem which was encountered involved the choice of which 
models to include in the set of possible models. The aircraft, of course, 
responds in a continuous way to parameter changes, and model set selection 
was found to have a considerable effect on perf ormance . 

It is issues such as these which have motivated our study. Our 
goals have been to understand the characteristics of the MMAC method 
and develop a useful design methodology based on the MMAC algorithm. 

1.3 Contributions of This Thesis 

This thesis presents the results of a detailed study of the MMAC 
algorithm. The major conclusions of this study, which are detailed in 
Chapter 6, are of two basic types. First of all, there are specific 
conclusions which the analysis of Chapters 4 and 5 yield. However, 
since the analysis of these chapters relies heavily on a special case, 
the results are of relatively little direct applicability. However, 
they are indicative of the types of and basic causes for behavior which 
has been observed in more general situations. Thus, of possibly greater 
importance are the qualitative conclusions which result when the specific 
conclusions are extrapolated to general problem structures. These 
qualitative results are also detailed in Chapter 6. 

The following are the major results of this thesis. 

1. The neutral stability of the MMAC method is established. 

2. Conditions which guarantee state convergence are derived. 

3. The MMAC algorithm consists of a bank of Kalman Filters, 
each corresponding to a hypothesized linear time- invariant model of 
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the system generating the observations (referred to as the true system) . 
The outputs of each filter are- used to calculate the aposteriori prob- 
ability that each filter “matches the true system- The’ feedback control 
is then given -by the probabilisticly weighted average of the controls 
calculated for each hypothesized model. (See Chapter 2 for a complete 
discussion) . This results in a highly non-linear closed loop system. 
However, if the probability is constrained to be constant , then the 
closed loop system becomes linear time-invariant. In this thesis it 
is shown that even if the closed loop system is unstable for all 
constant values of the probability, the overall system may have a 
bounded (in fact "hyperbolicly stable" - see Section 4.6) response. 

t 

For a special case, specific conditions are derived. 

4. The effects of numerical roundoff are examined and a form of 
implementation which behaves well is proposed. 

5. The specific results of this thesis can be used to predict the 
qualitative behavior of MMAC systems which have a more complex structure 
than the ones analyzed herein. 

4.1 Overview of this Thesis 

In Chapter 2, the MMAC algorithm is introduced in the form which 
will be used in the reaminder of this thesis. Both the continuous and 
discrete time versions of the algorithm are presented, although the 
discrete version is employed for the majority of the analysis. 

In Chapter 3, the canonical problem which forms the focus for the 
entire thesis is introduced. This problem, about which various structural 
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assumptions are made, has been carefully selected m order to provide a 
situation which retains many of the qualitative properties (such as the 
existence of oscillations) which have been observed m more general 
problems while simultaneously being amenable to detailed analysis. The 
remainder of the chapter contains a discussion of the various types of 
behavior which have been observed in simulations of the canonical problem 
along with sample simulation results to illustrate each behavior. The 
basic responses are shown to be of three types. In the first, termed 
"exponential" or "geometric", the states are geometrically stable. In 
this case, all ox the states are decreasing for all constant values 
of the probability. The second type termed "oscillatory" results in the 
probability oscillating between zero and one, which m turn results in 
the states exhibiting an oscillatory behavior, alternately increasing 
and decreasing along with the probability. The third type termed 
"mixed", results m a behavior which, depending on the magnitude of the 
initial conditions, exhibits either an oscillatory or exponential 
behavior. These simulations have been used to motivate the analysis 
in the remainder of the thesis. 

Chapter 4 contains the majority of the analysis of the MMAC method. 
In this chapter, each type of behavior described in Chapter 3 is analyzed 
in order to yield an understanding of the underlying causes of the 
behavior. This results in conditions which guarantee the existance of 
the exponential mode. Furthermore, various approximations are used to 
characterize the mag or modes of behavior which lead to conditions for 
the presence and absence of each behavior for the special case of the 
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canonical problem. Beyond this, more general qualitative results are 
also presented. 

In Chapter 5, various aspects of the problem of implementing the 
MMAC method using a digital computer are discussed. By far the most 
important of these is the modification to the analysis of the oscilla- 
tory behavior when the finite precision nature of the computer becomes 
a factor in the behavior. Also discussed are the effects of using 
each of the various forms of the equations as far as numerical accuracy 
is concerned. The chapter concludes with the proposal of a new form 
of the equations which is believed to allow the designer greater latitude 
in design without encountering numerical problems. 

Chapter 6 contains a discussion of various ad hoc design modifica- 
tions which have been proposed in order to overcome the shortcomings of 
the MMAC method. Also included there is a summary of the conclusions 
of this thesis as well as some suggestions for future research. 

1 . 5 Notation 

The following is a brief list of the notation employed in the 
thesis. Except for the notation for the components of the residual 
of a Kalman Filter, all are believed to be standard. 

Matrices are represented' by upper case letters which are 
underlined. 

Vectors are represented by lower case letters which are 
underlined. 

Scalars are represented by lower case letters ( not under- 
lined) . 
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A' 


I I A|| 

X(A) 
N (A) 




Ax(k) 


the transpose of A. 

the norm of a matrix = Max A. (A 1 A) 

the set of eigenvalues of a matrix. 

the Null Space of a matrix. 

the 0 -norm of a vector x given by 

/x'0 ^x for 0 ^ positive definite. 

the norm of a vector given by Vx'x 

magnitude of a real or complex a. 

change in x(*) from k-1 to k. 


N 

H - Product for i=l to N 

i=l 


N 

s 

i=l 

x 


x. 

i 



i 


(a,b) 

[a,b] 


- Sum for i = 1 to N 

- vector of true state variables 

- scalar true state component i 

- scalar j — component of r 

-> —i 

- open interval between a and b. 

- closed interval between a and b. 


#§£■35355 



CHAPTER 2 


REVIEW OP THE MMAC METHOD 


The purpose of the present chapter is to introduce the Multiple 
Model Adaptive Control (MMAC) algorithm. A full discussion will not be 
given as that is available from other sources [7, 12]. 

The MMAC algorithm is composed of two parts. The first, which 
performs an estimation/identifxcation function , is similar to a Maximum 
Aposteriori Probability (MAP) algorithm [15] which is discussed in Section 

2.1 for the discrete time case. The MAP algorithm is structured as a 
bank of Kalman Filters with some decision logic. The second part, which 
is cascaded with the MAP- like algorithm, is a control computation which 
is discussed in Sections 2.2 and 2.3 for the discrete time case. 

The remaining sections of this chapter contain a development 
and discussion of the special forms of the equations for the MMAC algorithm 
which prove useful in the remaining chapters of this thesis . 

2.1 Maximum Aposteriori Probability (MAP) Identification 

2.1.1 The Kalman Filter (KF) 

Assume that a linear, time-invariant (LTI) discrete time system is 
given by: 

x(k+l) = A sc(k) + B_u.-(k) + C_(k) (2.1a) 

with observations : 

y(k) = Cx (k) + n(k) (2.1b) 
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where x(k) is the n-dimensional state vector, u(k) is the m-dimensional 
control vector and y_(k) is the p-dimensional observation vector. The 
noise sources j^(k) and 7]_(k) are taken to be zero mean white Gaussian 
noises of covariances = and ^respectively. The matrices A (nxn), 

(nxm) , and <pxn) are the system, input and output matrices respectively. 
We will use the notation 


(A, B, C) (2.2) 

to refer to the above system. The system (A, B, C) which generates the 
observations ;g_(k) will be called the time system . 

In practice , the actual values of the matrices A, B , £, £ and ^ are 
unknown. However, estimates of these parameters are often available from 
a knowledge of the system. Thus, (A^, B^, C^) will be used to denote 
the i model of the system, given by: 


x(k+l)' =s A.x(k) +■ B.u(k) + £(k) (2.3) 

— -a— —i — 

y (k) = Clx(k) + n_(k) . 

For the purposes of the present study the values of ip. and will not 
vary from model to model, although extensions to that case could be con- 
sidered. 

It is well known that the steady state Kalman Filter (KF) [19, 20] 
which estimates the state x(k) , based on the model 



is given by 



- 21 - 


ORIGINAL PAGE IB 
OF POOR QUALITY 


x. (k+1) — A.x. (k) + B.u(k) + H.r. (k+1) 


r(k+l) = y(k+l) - C [A.x. (k) + B u (k ) ] 


, . th 

where H. is the Kalman gain for the i — model: 


(2.4) 


H. = S. CMp 


-1 


(2.5) 


and Z. is the nxn solution to the familar Riccati equation 
~i 

-1 ’ -1 
2. =* [c'Tp c. + (a.e.a : + =) ] 


( 2 . 6 ) 


The state estimate x. (•) is n-dimensional, while the residual vector. 


—l 


r. (-) , is p-dimensional. 


2.1.2 Properties of the Kalman Filter 

The Kalman Filter has many interesting properties , (see for example 

[19]) , a few of which are useful in understanding the MMAC method. These 

are now discussed. A few definitions are useful. 

tii , 

Definition 1: The i Kalman Filter is said to be matched if the matrices 

used in the filter design (i.e. , the model) and the manrices of the true 

system are identical; that is , if A = A B = B_ and C. = C. 

1 r “ f 1 * '1 

Definition 2 : The Filter is called mismatched if it is not matched. 

til -i 

Property 1 : If the i Kalman Filter is matched then in steady state: 


E{r. (k) } = 0 


E jr\ (k)9 i 1 r^( e )|= l6(k-e) 


(2.7a) 


(2.7 b) 
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where <S C - ) is the impulse function, e{*} denotes expectation, and 8^ is 
given by ‘ 

0^ = jjj. + C_£ ^C' . (2.7c) 

Furthermore, if the i^ 1 filter is matched, then x^(0 gives the optimal 
estimate of the state x ( * ) . 

til 

Property 2 ; If the i Kalman Filter is matched, then the probability 
density function for the residual ( ♦ ) is given by [20]: 

- ~ [rieT 1 !] 

p(r.) = (3.e ^ (2.8a) 

— i i 


with 


( / — ™ — rr 1 


s i = H 2 } lii 1 / 

tli 

Property 3: If the i Kalman Filter is mismatched then 

(2.8b) 

E{r. (k) } - r (k) 
— i -i 

(2.9a) 

e{ ( r. (k) - r. (k))07 1 (r U) - x (£))'}» F. (k-&) 

—i —i ~i —i - i — i 

(2.9b) 

where r in general is nonzero and 0^ is a function of the 
noise covariances (see [28]) with, in general, 0^(0) > I_. 

system and the 

It follows that if two filters (one matched and one mismatched) are 
computed then 

E{r^ (k) 0. (k) } < E{r_^ (k ) ^—2 ^ ^ 

(2.10) 


where r is from the matched filter and r_ is from a mismatched filter 

a. ^ 
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[16]. Thus, if equation (2.8) is evaluated for both a matched (1) and 
a mismatched (2) KF, then 

s{p (r^) } > E{p(r 2 )> . (2.11) 

2.1.3 The Probability Equation 

Assume that we are interested in identifying a dynamic system from 
its outputs and that the true system is known only to be one of N speci- 
fied models. Baye’s rule and Properties 1-3 imply [15, 7] that the 

iili 

probability that the l model matches the unknown system (P^) is given 
by 

P. (k) p (r (k+1)) 

P.(k+1)=*-— — (2.12) 

1 N 

E P . (k)p (r (k+D) 

j»l 3 ^ 

where p (r_^ (k+1) ) is given by equation (2.8). The struct lire of the re- 
sulting system is shown in Figure 2.1. 

2.2 Adaptive Control 

If one knew wiuh certainty which model matched the true system, it 
would be a simple matter to design a controller using any of the standard 
synthesis techniques. Therefore, one reasonable way to determine a 
control law for the unknown system is to probabilisticly weight the con- 
trols which would be used if one assumed that one of the models was 
correct. That is, let 

N 

u(k) = E P (k)u. (k) 

— . , i — i 

i=l 


(2.13) 
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where p 00 is given by equation (2.121 and u^OO is the control which 
one would apply if model i were assumed to be matched to the true system. 
It will be assumed that 

u^ = GjsL (2 .14) 

although this is not necessary and is further discussed in Section 2.3. 
Figure 2.2 thus summarizes the Multiple Model Adaptive Control method. 


2.3 The Control Law 

The MMAC method as developed in the literature [7, 16] has assumed 
that the linear quadratic (LQ) methodology is used to design the controller. 
Thus, the feedback gains of equation (2.14) are chosen to minimize 

CO 

J (u) = 2 (x' (k)0 3c (k) + u' (k)R;u (k) ) . (2.15) 

k=l 

It can be shown [21] that the optimal gains are given by 

% “ + (2.16) 


where K. is the solution to the steady state Riccaui equation 


K ± = Cn + - A] Kj. 3 ^ + R] 


Although to date, all references to MMAC have made use of the control 
law (2.16), this is not a necessary part of the method. Thus, any control 
law which gives good results for the respective model may be used. How- 
ever, there is a strong interaction between control law choice and 
adaptive performance due to the feedback, about which very little is known. 




Pig. 2.2 Adaptive Control Structure 
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For the purposes of this study we will for the most part restrict our 
research to control laws based on the linear quadratic (Equation 2.15) 
form. 

2.4 Comments on MMAC 

A few comments are in order. 

1) The MMAC algorithm, as shown by Willner [7] is suboptimal even 
for the problem originally posed. Willner was able to show that the 
algorithm is optimal for the final step of the dynamic programming 
algorithm [7] but was unable to continue the calculation backward in 
time. 

2) As posed above, the true model is assumed to belong to a finite 
set of known models. Lainiotis 19, 13] has discussed the infinite set 
case but concludes that a finite approximation is then required for use 
in applications . Thus , for most real problems when the true model may- 
take values from an infinite set, a further suboptimal approximation is 
required. The results of Baram 115] (which apply only to the open-loop - 
case) may help in discretizing the model set. 

2 . 5 Continuous Time MMAC 

As developed in the preceeding sections, the MMAC method has been 
based on discrete time systems. For analysis purposes, it will be useful 
to consider the related continuous time problem, largely because of the 
simplified form of the probability equation. 

The complete equations for the method will not be given here as they 
are the continuous time regulator and continuous time Kalman Filter 
equations [19, 6] with a set of equations for the probability of each model. 



Dunn [22] has shown that if P (t) is the probability of the i— 
model at time t then 


N N 

P. (t) = P.Cc.x. - Z P C.x.)0 _1 (y - E p.C.x ) (2.18) 

1 1 j=l - - D=1 3-3-D 

where Ch , x^ and y are defined analogously to the discrete time case and 
0_ is the observation noise covariance. The property which makes Equation 
(2.18) useful for analysis is the absence of either exponentials or a 
denominator. These equations are examined further in the two model, two 
state case in the following section. 


2.6 A Special Case 

Of special interest to the research at hand is the case when N=2 
(i.e., the two model case) with full state observation. Furthermore, we 
shall assume that the input (3) matrix is the identity. 

The special case captures all essential features of the problem that 
we are interested in examining without adding unnecessary complexity, and 
it forms the central focus for this research. 

2.6.1 Discrete Time Case 

In the discrete time case, the true system is then given by 


x(k+l) = Ax (k) + u(k) + j£(k) 
y (k) = x(k) + T)(k) 


(2.19a) 


with models 


Model 1: £ (k+1) = A^tk) + u(k) + £(k) 


(2.19b) 



(2.19c) 


Model 2 : x_ 2 (k+1) = & 2 ± 2 Otl + u(k) + £Cfc) 

y_(k) = x 2 (k) + n_(k) 

Thus, the MMAC method reduces to the following set of discrete time 
equations 


x(k+l) = Ax (k) + u(k) + (k) 

Y (k) = x(k) + n(k) 

x, (k+1) = A^x^Ck) + u(k) + (k+1) 

r_ (k+1) = y_(k+l) - A (k) - u_(k) 

x_ 2 (k+1) = —2—2 ^) + — + -2~2 (IC+1J (2.20) 

r 2 (k+1) = y_(k+l) - A Xg (k) - u(k) 

u(k) = -P^x^k) - ^ 2 — 2—2 ^k) 

P (k) p (r ) 

P (k+1} = — i 

1 P 1 (k)p(r 1 ) + P 2 (k)o(r 2 ) 

- | [r.'e^r.] 

p(r x ) - B ± a 2 ^“ 1 “ 1 . 


The major goal of the research has been to understand such qualitative 
properties of the MMAC method as stability. Since the phenomenon which 
have been observed are largely due to the nonlinearities of the system 
and occur even in noise-free simulations , the noise terms t, ( • ) and n ( • ) 
add an unnecessary complexity. Thus, these terms are ignored in the 
remainder of this thesis. It should be noted that the KF's, designed 
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assvtming the noises to be present, are an integral part of the MMAC method 
and are therefore retained. Thus, the principle focus of the work is to 
examine the deterministic properties of Equation (2.20) and the correspond- 
ing continuous time system discussed in the next subsection. Since the 
sum of the probabilities of the models must always be 1 , it is known that 

P 2 (k) = (l-P ]L (k)) Vk. 

Then, rewriting Equation (2.20) in terms of the residual r. (k) we see 

1 

that the equations of the MMAC method can be summarized by 


x(k+l) = Ax (k) + u(k) 


r^ (k+1) = (A-A 1 )x(k) + A^ClyH^Jr^Ck) 


£ 2 (k+1) = (A-A 2 )x(k) + ^ 2 Cly_H 2 )r 2 (k) 


( 2 . 21 ) 


u(k) = - <p G + (l-P x )G 2 )x(k) + P X G (l-H^Jr + (1-P 1 )G 2 (I-H 2 )r 2 


P x (k) p(r x ) 


p i ( k+!) p {k) + ( i- Pi ( k }) p( r 2 ) 


p(r ) = B.e 

—a. i 


~ [r:e Vi 

2 —i— -a 


It will be useful notationally to combine the state and residual 
equations into one vector equation. Therefore, we define the vector 


w* (k) = [x‘ (k) , £ x (k) , £ 2 (k) ] ' 


Thus, Equation (2.21) can be rewritten as 



ORIGINAL PAG® IS 
OF POOR QUALITY 
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or 


w(k+l) 


A-P,G, - (1-P ) G_ 

P, G, (.1— H ) 

(1-P,)G„ (I-H„) 

— 1—1 1 —2 

1—1 -1 

1 —2 2 



0 

(A— A ) 

0 

A_ (I— H ) 

L — 2 


—2 2 _! 


(2.22b) 

w(k) 


w(k+l) = A(P^)w(k) 


with 


p (k+1) 


p (k)p(r ) 

P-^k) p(r 1 ) + (l-P^Ck) )p(£ 2 ) 


(2.22b) 


These equations , along with their continuous time counterparts , form the 
basis for the research which has been undertaken. 


2.6.2 Continuous Time Case 

Similar assumptions to those previously presented for the discrete 
time case can be made in the continuous time case. Here, we will again 
restrict our attention to the two model, state feedback case with B=I. 
Thus, the MMAC method reduces to the following equations 

x(t) - A 31 (t) + u_ (t) + £(t) 

ytt) = rCt) + Tjtt) 

x^Ct) = A * (t) + uCt) + H^r^Ct) (2.23) 

r^Ct) = y_(t) - (t) 

x 2 (t) = A 2 £ 2 (t) 4 - uCt) + H 2 £ 2 ^ 
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r 2 (t) = z (t) - x 2 (t) 


u{-t)- = -Pj-G^Ct) - P 2 £ 2 i 2 (t-> 

p x ct) = p 1 ct) p 2 ct) ig ± ct) -x 2 (t) j 1 e " 1 ix-p i ct)S r ci-p 1 (t) )g 2 ct) j 


with P 2 (t) = 1-P 1 (t). 

For the same reasons as presented m the previous section the noise 
sources £_(t) and jn(t) will again be ignored for the majority of the 
proposed research. It should be noted that under this condition, the 
residual r^ (t) exactly equals the estimation error x(t) - x^(t). Thus, 
equations (2.23) can be rewritten as 

V 

x(t) = Ax (t) + u(t) 


r (t) = (A -H^r^t) + (A-A^xO:} 62.24} 

r 2 (t) = (A 2 -H 2 )r 2 (t) + (A-A 2 ) x(t) 
u(t) = -P^tjG^xtt) - <t> ) - (1-P (t)) G 2 (x(t) - £ 2 Ct)) 

P 1 (t) = P 1 (t) [1-P^^ (t) ] [r 2 (t) - r^tt) 3 ' £ i lP- L (t)r 1 (t) + (1-P 1 (t) )r 2 (t) ] 
Thus , if we again combine equations by defining the variable 


x(t) 


w(t) = 


r^tt) 

r 2 Ct) 


Equation (2.24) can be rewritten as 
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w(t) = 


a - e A - u-p^jj j a - Ei )g 2 




(A-A 2 ) 


t 


v n. 

- -2 


A„ - Ih 

-2 


w(fc) (2.25a) 


or 


w(t) = SCP^lwCt) 


with. 


3?l (t) = P 1 (t) [1-P 1 (t)]]r 2 (t) 


(2.25b) 

(t) ] 07 1 lP 1 (t)r 1 (t) + (l-P 1 Ct) )r 2 (t) ] 


2.7 & Change of Variable 

A change of variable which proves useful in the sequel is to let 


q = 2P - 1 . (2.26) 

Making this change in Equation (2.22) yields 

w(k+l) = A(q) w(k) (2.27) 

( 1+q) p (r^) - (1-q) p (r ? ) 
q(K+1) = “(l+q)p (£")“+ (l-q)p(r 2 ) 


where 


A - 


A(q) = 


ia4,)G 1 -i( 1 -<l)£ 2 


A - 


L 


h.- b . 2 


yCl+q) G CXj-H^) 


2 2 ^ I_ — 2^ 

0 

V!-h 2 ) 


In continuous time the corresponding change to Equation (2.25) yields 
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w(t) = A(q) w (t) 


q(t) = i(l-q 2 ) rr 2 -r 1 ]'0" 1 [(l+q)r ;L + (l-q)r 2 3 



A - yd+q) - jCl-q) £ 2 


|-a-q)£ 2 

A(q) = 

A - A 

h'h 

£ 


A — A. 

Q 

A„ ~ H- 


~ —2 


—2 — 2 1 


(2.28) 


Note that the same basic notation is used for the A matrices of Equations 
(2.22), (2.25), (2.27) and (.2.28). However, no confusion should arise 
as the meaning should always be clear from context. 


,2 . 8 A Useful Definition 

A further concept which will prove useful can be seen by considering 
either the continuous or discrete time problems as summarized by either 
Equation (2.25) or (2.22) 


Continuous time: w(t) = A(P )w (t) 


Discrete time: w(k+l) = A(P^)w(k) 


(2.25). 

( 2 . 22 ) 


where 


Continous time : 


A(P 1 ) = 


Discrete time ; 


- _P l- 1 *' Cl " P 1 ) -2 


— ~ — i 

L - “ -2 

A-^Gl-Cl-Pl)^ 




A - 


P A 


-i'5l \ ° 




-1 






(1 -VS2«-S^ 




4 
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Thus, it is clear that if P^(-) is held constant, then Equations (2.25) 
and (2.22) describe linear, time- invari ant dynamical systems in the 
variable w' = [x' , r£, r^3 ' . Thus, Equations (2.25) and (2.22) with 
P constant will be referred to as the linear, time invariant system for 
fixed p , or LTI for fixed P. As will be seen in Chapters 4 and 5, many 
of the major properties of the MMAC method can be expressed in terms of 
the properties of the LTI system for fixed P as a function of P. 

A second important point to note is that in applications and 
simulations of the MMAC system, another parameter becomes important due 
to the finite precision usually available for computation. Prom the 
equations of this chapter, it is clear that the MMAC method has the 
property that 

P. (k) =0 =?>P. (•) =0 for all future times. 

i i 

Such a situation is usually to be avoided; - as it reduces the flexibility 
of the algorithm m a number of applications. For example, the parameter 
of the true model are often changing in adaptive control situations , and 
one would like to require P to be nonzero for all time so that the MMAC 
algorithm can respond to such changes. Thus, one almost always applies 
an additional constraint on the probability such as 

P. (k) > P, fi, . 

i — lim 

-50 

This has been done throughout this study with a value of P - 
The effect of such a limit is examined m detail in Section 5.1. 




CHAPTER 3 


QUALITATIVE RESULTS 


In order to guide and motivate the research, we have examined a 
prohlem consisting of a system with two independent states and two models. 
This system, while simple, captures many of the basic issues which are 
important to the method and sheds light on the fundamental problems in- 
volved with the MMAC design. The problem is formulated in the next 
section of this chapter. The remaining sections contain simulations 
which demonstrate the various types of behavior which have been observed. 
A discussion of the important properties is included. 

3.1 Problem Formulation 

In most applications of an adaptive control algorithm, a detailed 
analysis of the behavior of the algorithm has proved intractable. This 
is especially true for the MMAC algorithm. However, m the case of 
MMAC, it is possible to find a simple exaitple problem that lends itself 
to analysis and simultaneously maintains the basic properties exhibited 
in simulations of the more general systems. Thus, for the purposes of 
this thesis, a sample problem structure has been chosen which displays 
what we feel are the critical characteristics of the method and which 
is also amenable to detailed analysis. 

The chosen true system to be .controlled is given by: 


x(k+l) 


a 

0 


0 

a 


x(k) + u(k) + £(k) 


y_(k) = x(k) + rKk) 


(3.1) 


- 36 - 
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or a continuous time system of the same structure and A matrix 

x(t)' + u(t) + ?(t) (3.2) 

y.(t) = x(t) + n_(t) 

where a takes on values xn the range [0, 2]. The discrete time system 
is useful for simulation studxes and most of the analysis. However, use 
of the continuous time system provides greater xnsight for the lineariza- 
txon results in Section 4.1. The set of models for either the discrete 
or continuous time case is gxven by 



Model 1: (A^, I, I) 

. Model 2; (A 2 , jl, I) 

wxth 

■ /v » *1 

a 0 
0 a_ 



where a takes on various values from 0 to 1.5. The parameters a and a 
are varied to obtain different responses from the overall system. 

For the purposes of KF design, the noise sources _?(*) and ]]_(•) are 
assumed to be zero mean, white and Gaussian wxth covariances of I<5(*) 
where I_ xs the 2x2 identity matrix. This structure results in diagonal 
Kalman Gaxn matrices such that 



'h 0' 


A 

"h O' 

H, - 


= 


-1 

A 

—2 



.0 h_ 


.0 h_ 
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where h is the gain corresponding to a and h is the gain corresponding 
to a. Further, the control law weight Q. and R. are chosen diagonal such 

”1 ■*” X 

that 



0 

1 


1 — 

o 

<tr* 



I! 

C}J 

Oil 



_ /\ 

|_0 q_ 


0 q. 


This results in a diagonal control gain matrix 



The structure described above will be referred to as the canonical problem. 

The emphasis of this study has been on the examination of such quali- 
tative properties as stability . Thus , we ignore the noise sources both 
in the simulation and in the analysis of the properties of the MMAC 
method. That is , in both the analysis and simulations , £ ( * ) and p ( ■ ) are 
set equal to zero. Note, however, that the KF's designed with the assumed 
non-zero noise sources are retained. It is clear that noise can have 
a major effect on a system [29]. However, as we will see, many of the 
properties of the MMAC method are due to the nonlinearities of the 
probability equation. Thus, it is our feeling that an anlysis of this 
noise-free case is of considerable importance. 

A few comments are in order concerning the the choice of the example 
problem. First of all, the system is chosen to be the simplest possible 


and still capture the important phenomena. Thus, we selected a two 
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state, two model system. The choice of model is, admittedly, somewhat 
extreme in the degree of symmetry and mismatch between the models and 
the true system. However, this choice has been deliberately made in order 
to investigate phenomena which have been observed in actual applications 
[23J . Thus, it is felt that the analysis of this problem has yielded 
insight into the more general case. 

3.2 Basic Responses 

Various types of responses have been observed in simulations of the 

system presented in Section 3.1. Table 3.1 details the parameters used 

for each simulation. The remainder of this chapter presents examples 

of each of the major modes along with a discussion of the important 

characteristics of each. The simulations, which have been performed on 

an IBM 370 computer using double precision FORTRAN, have been used to 

motivate and guide the analysis of Chapters 4 and 5. Xn fact, the types 

of analysis described in these chapters and, in particular in Section 

% 

4.6, were to a large degree suggested by the simulations discussed m 
this chapter. 

For each simulation presented, three plots have been included. The 
first is a plot of the probability of model one (P^) versus time. The 
second is a plot of the two true states and versus time. The third 
is a plot of the quantity ln^x^. This quantity has been found to be 
indicative of the stability of the closed loop, nonlinear system. This 
variable is further discussed m Section 4 . 6 where it is linked with 
the concept of "hyperbolic stability”. 
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Case # 

a 

a 

q 

g 

A 

q 

A 

g 

h 

fi 

Figures 

la 

2 

0 

1-02 

1.62 

i 

0 

.809 

.5 

3.2 

lb 

2 

0 

0 

1.5 

l 

0 

.809 

.5 

3'. 3 

lc 

2 

0 

* 

1.4 

l 

0 

.809 

.5 

3.4 

2 

1.5 

1.0 

1 

1.09 

i 

.618 

.7245 

.62 

3.5,6 

3 

.9 

0 

1 

.538 

i 

0 

.597 

.5 

3.1 


TABEE 3.1 

Parameters of Sample Cases 


*This value of the control can not result from an LQ design 
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3.2.1 Exponential Mode 

The first type of behavior is one in which the states of the true 
system increase or decrease in an exponential manner (see Figure 3.1). 

This type of behavior appears to arise in two situations. First, 
examination of the equations of Chapter 2 indicate that if the KF resi- 
duals and r^ are equal with the true state components equal to each 
other, then P^t) = 0 or equivalently, P^Ck+l) = P^k). Thus, if the 
system is symmetrically initialized (that is, the two true states are 
equal with r ^ = _0 and P^tO) = .5) then (t ) = .5 ¥ t or equivalently 

P^ (k) = . 5 ¥ k . The closed loop system is then time invariant and sta- 
bility analysis for the resulting LTI system follows as usual. Note that 
the resulting system can be exponentially stable, neutrally stable or 
exponentially unstable depending on the control gains and G^* al- 
though this is clearly a singular condition, it nonetheless is important 
from an applications point of view because one commonly attempts to 
initialize the system with equal probability and with r_ = r^ (i.e., 
x = x ). An analysis of this mode is presented m Section 4.1. 

A somewhat similar type of behavior occurs when the LTI system for 
fixed P is stable for all P. This, of course, is a fairly restrictive 
condition which m essence requires extreme robustness of each controller. 
However, such a condition does imply exponential stability. This occurs in- 
dependently of how the probability behaves. Note that this is a non- 
trivial result since A(P(k)) is time-varying because P(k) is. Section 
4.2 contains the analysis of this situation and Figure 3.1 is an example 
of the type of simulation results obtained. 
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3.2.2 Oscillatory Reponse 

Probably the most unusual behavior which has been observed both in 
the current work and' in applications of the MMAC algorithm [23] has been 
an oscillatory response in which the probability jumps between near- 
zero and near-one* in what appears to be (but strictly speaking is not) 
a limit cycle. 

Figures 3.2 through 3.4 are examples of this type of behavior. Figure 
3.2 is a case m which ’ the peaks of the state traj ectories are approxi- 
mately constant while the period of oscillation increases. Figure 3.3 
is a case in which both the peaks and the period are constant while in 
Figure 3.4, the system is unstable . These three cases are obtained by 
changing the value of the control gain. Note that each gain would yield 
stable behavior for the model used m its design. 

It is interesting to note that the states of the system are also 
highly oscillatory, it might be expected that the plant dynamics would 
smooth the rapid probability transitions to form a smoother "average” 
state response. However, as shown in Chapter 4, the oscillatory state 
behavior is a direct consequence of the model mismatch problem m the 
MMAC algorithm. 

The reasons for this oscillatory behavior, which are discussed in 
Sections 4.4 through 4.6 and again in Section 5.1, are closely related 
to the fact that neither of the hypothesized models individually yields 

*This is the one set of simulations in which the lower limit on the 
probability (see Section 2.8) of 10“ 50 is achieved. Section 5.1 dis- 
cusses the analysis in this case in detail. 
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a control which stabilizes the true system. Basically, the adaptive 
control law then attempts to achieve stability by alternately controlling 
each mode of the system. This results in unusual behavior which may be 
bounded, (figure 3_._2 and 3. 3 ) or unstable (Figure 3.4). Chapters 4 and 
5 contains analysis which gives conditions for each type of behavior. 

The principle conclusion is that stability will result when the controller 
at any time stabilizes some modes faster than it destabilizes the re- 
mainder. Thus, for the two state problem, a state must be reduced more 
when the controller stabilizes it than it is increased when the control 
results in unstable behavior for that state. 

One variable which appears useful in characterizing the above condi- 
tion is the product of the true states , that is , the quantity 
variable is plotted for each simulation (i.e., see Figure 3.2(c)). a 
complete analysis of the properties of this variable is given in Section 
4.5. While not strictly a Lyapunov function, examination of this variable 
provides a type of analysis which permits a characterization of the sta- 
bility behavior of the MMAC method. Furthermore, it captures the important 
observed characteristics of the simulations. For example, the analysis 
of Section 4.6 uses this variable to predict the stability of the three 
cases of Figures 3.2 - 3.4. The connection with stability can best be 
seen by considering a plot of versus x^- Figure 3.5(a) is such a 
plot corresponding to the simulation of Figure 3.2. The state trajectory 
for this example tends to look like a family of hyperbolas 

bt 

x = e 
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If the overall system is bounded, the trajectory then changes as in 
Figure 3.5(a). Likewise, unbounded behavior occurs whenever b is 
greater than zero (see Figure 3.5(b)). 

It should be noted that simply bounding the product of the states 
is not sufficient to guarantee the boundedness of the individual states 
since the states could still go to infinity along a hyperbola. However, 
an analysis of the probability equation (Sections 4.4 and 5.1) provides 
bounds on the values of the state at which probability transitions will 
occur. This leads to the conclusion that the peaks of the state tra- 
jectories will be bounded whenever the product of the states is bounded. 
Thus when a constant value of P results in a stabilizing control, the 
system may still be bounded but not asymtotically stable. Further, 
when the system is bounded and b is less than zero, then the period of 
oscillation is increasing since the stable mode decreases more than the 
unstable mode increases and the peaks of the curves are constant re- 
sulting in the unstable mode having to increase more on each cycle of 
the probability. 

3.2.3 Mixed Case 

When the LTI system for fixed P is stable for some P but not for 
all P, depending upon the initial conditions, we can obtain simulations 
which exhibit characteristics of either of the preceding types of be- 
havior m that either an oscillatory or an exponential response may be 
observed depending on the initial conditions. For large initial con- 
ditions the response will initially be oscillatory but usually will 
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finally decay to a constant value of P . * Furthermore , this limiting 
value of P will be such that A(P) is a stable matrix. Figure 3.6 is a 
simulation of the Case 2 configuration with large initial conditions'. 

When the initial conditions are small, a second type of response 
occurs. Assume that A(P) is a stable matrix for P £ [l/2-£, 1/2 +£ 3 
and unstable otherwise. If P<0) = 1/2, then there exists some non-zero 
w(0) such that j jw(k) | [ decays sufficiently quickly so that the resulting 
change in P does not. take it beyond 1/2 + £. This is a direct consequence 
of the fact derived in Section 4.1 that the change in the probability is 
proportional to | [wj | while the change in w is proportional only to 
| |w| | . Figure 3.7 is a simulation of the Case 2 configuration for small 
initial conditions. Note that the probability merely makes a small 
jump and that it shows little tendency to return to 1/2. Section 4.7 
details a procedure which can be used to estimate the range of w(Q) 
which results in this non-oscillatory behavior. 

3.3 Summary 

The preceedmg sections have given an overview of the types of 
behavior which the MMAC method can produce. Table 3.2 summarizes the 
major characteristics of each. Each type of behavior requires a dif- 
ferent form of analysis in order to understand the dominant . effects . 

This analysis is given in Chapters 4 and 5. 

Exponential behavior occurs primarily when the basic characteristics 
of the closed- loop system are independent of the probability, either 

*Whether or not the oscillations always die out in this case remains 
an open question. 
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Values of 
3? yielding 


Section 

Fig. 

Stable LTI 
System 

b. 

i 

bounded 

■p 

period 

of oscillation 

3.2.1 

3.1 

All 

decreasing 

Yes 

NA 

3.2.3 

3.5 

3.6 

Some 

decreasing 

Yes 

increasing 

3.2.2 

3.4 

None 

increasing 

No 

•k 

3.2.2 

3.3 

None 

constant 

Yes 

constant 

3.2.2 

3.2 

None 

decreasing 

Yes 

increasing 


N A = resulting behavior not periodic 
* = period approaches a constant 


TABLE 3.2 


Principle Modes of Response 
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because the probability is held constant by the dynamics or because all 
states either increase or decrease regardless of the value of P. Sections 
4.1. through 4.3 provide significant insight in this case. Section 4.1 
contains the analysis of - .the local behavior, concluding that the prob- 
ability equation is neutrally stable because the rate of change of the 
probability is proportional to the square of the system states. Section 
4.2 considers the only case in which exponentially stable behavior can 
be guaranteed. This is shown to occur when A(P) is a stable matrix for 
all values of the probability. Section 4.3 provides an asymtotic result 
when one of the models matches the true system. 

The oscillatory mode of behavior is analyzed m detail in Sections 
4.4 through 4.6 and includes criteria as to when the overall system will 
be stable or unstable. It is shown that the important parameter de- 
termining stability is the relationship between the stabilizing time 
constants of each state. Thus, if the controller stabilizes one state 
faster than it destabilizes the other, the overall system will be stable. 
Section 4.4 contains a detailed analysis of the probability equation 
which shows why state oscillations are an integral part of the oscillatory 
behavior. Finally, Section 4.7 discusses the reasons for the mixed type 


- of behavior. It concludes with a procedure for determining the range of 
initial conditions such that non-oscillatory behavior will occur. 



CHAPTER 4 


STABILITY ANALYSIS 


It is the purpose of the present chapter to provide the analysis 
necessary to understand the basic stability properties and qualitative 
behavior of the MMAC method- At present, there is no single stability 
result which totally describes these properties. Thus, it becomes 
necessary to examine several approaches, each of which adds no the over- 
all picture, but with none providing a whole view. Furthermore, it 
often is necessary to combine the results of differing methods to deduce 
a single property. An example of this is the combining of the Lyapunov 
results of Section 4.5 with the analysis of the behavior of the probability 
equation m Sections 4.4 and 4.5 to arrive at a stability result which 
neither type of analysis alone could provide - 

Chapter 3 contains simulations of carefully selected special cases 
of systems controlled by the MMAC algorithm. The purpose of the present 
chapter is to attempt to provide an understanding of the properties of 
the MMAC method by first noting various features from the simulations 
and then , guided by the simulations , attempting to understand the features 
by analyzing some special cases. This results in some specific conclusions 
for the special cases and, more importantly, it yields considerable in- 
sight into the qualitative behavior of the MMAC system in more general 
situations. The major conclusions of this chapter are: 

1) At best, the MMAC system is neutrally stable about an equilibrium 
point m that the probability has no tendency to return to its initial 
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value following a perturbation. This is not an unexpected property for an 
adaptive controller. 

2) If A(.P) is an unstable matrix for p = 1/2 then the special case 
of Section 3.1 results in the probability oscillating. 

3) This oscillatory behavior results in the states being either 
bounded or unstable. Specific conditions for two special cases are given. 
Qualitatively, the requirement for the boundedness of the states is that 
the modes of A(P) which are unstable for P = ltP=0) must be dominated 

by the stability of the same mode when P = 0 (P=l) in that they must grow 
slower for P = 1 (P=0) than they decay for P = 0 (P=l) . 

4) The oscillations observed in Section 3.2 may be bounded even if 
no constant value of P results in a stabilizing control. The controller 
then attempts to achieve stability by alternately controlling each mode 
of the system. This alternating of controls also occurs for large 
initial conditions when A CP) is a stable matrix for some but not all 
values of P. 

As seen in Section 3.2, two major types of behavior have been ob- 
served in simulations : oscillatory and exponential (jion-oscillatory) . 

The analysis of the present chapter is aimed at understanding each. 

Thus, Sections 4.1 through 4,3 deal largely with the exponential behavior 
while Sections 4.4 through 4.6 consider the case m which the probability 
oscillates. Section 4.7 considers the case in which either type of 
behavior can occur and contains a discussion of the initial conditions 
that lead to each type of behavior. A detailed overview of the chapter 


follows 
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Section 4.1 through. 4.3 contain preliminary material dealing with 
the exponential case. For example, the neutral stability of the prob- 
ability exhibited in the simulations of Figures 3.6 and 3.7 is derived 
in Section 4.1. The other major conclusion of rhe section is that 
changes m the probability are proportional to the norm of the state 
squared while changes in the states are only proportional to the state 
to the first power. Thus, for small values of the state, the state tends 
to change faster than the probability. This further explains the switch- 
ing behavior exhibited m the oscillatory responses in Section 3.2. In 
Section 3.2.1, the case in which both models result m stabilizing con- 
trols for all values of the probabilities has been simulated. It is 
shown that under some conditions the probability dynamics can effectively 
be ignored in determining stability. Section 4.2 contains the analysis 
of this case and proves that a(P) must be a stable matrix for all P in 
order to guarantee that the states will be exponentially stable for 
all initial conditions. 

Section 3.2.2 includes simulations in which the probability tends 
to look like a switch, taking on values near one or zero but seldom in 
between. Section 4.4 contains a detailed analysis of the probability 
equation which leads to an understanding as to why the probability ^umps 
so rapidly. Sections 4.5 and 4.6 continue this analysis to present two 
methods of analysis for ascertaining when the oscillatory behavior will 
be bounded or unstable. These sections, which result m basically the 
same criteria, take different points of view; Section 4.5 contains 
approximations to the general solution over several time intervals while 



-70- 


Section 4.6 employs a function which resembles a Lyapunov function. 

From Section 3.2.3, it is clear that the size of the initial con- 
ditions partly determines whether the probability will oscillate or not. 
In Section 4.7 the situation is analyzed in detail and a procedure de- 
veloped which allows for the determination of the limits of the initial 
conditions which result in the probability not oscillating. 

The majority of the analysis in this chapter has been performed m 
the discrete time domain. In many ways this is to be preferred as this 
is the form in which the method is most often implemented m practice. 
However, some of the analysis is done using the continuous time version 
because the results are simpler and lend themselves to interpretation 
more easily. No claim is made as to the complete equivalence of the 
two forms. In fact, one could expect some differences between them due 
to the different assumptions on the availability of sensor data (i.e., 
discretly versus continuously available) . However, the qualitative 
conclusions of one are applicable to the other, as can be seen by the 
examination of the two sets of equations. It should also be bom m 
mind that the analysis of this chapter is aimed at understanding the 
phenomena exhibited in the simulations of Chapter 3. Thus, unless 
stated, it can be assumed that the canonical structure of Section 3.1 
is under consideration. However, the structure of the canonical prob- 
lem has been carefully selected to accentuate certain types of behavior 
observed m more general MMAC simulations [16, 23], Thus, the intuition 
and qualitative results are believed to provide significant insight in- 
to the behavior for more general systems. 
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It should be noted that throughout this thesis, the terms "oscillatory”, 
"non-oscillatory" and "exponential" refer to qualitative behavior of the 
probability . Although not observed m the simulations in Chapter 3, it 
is possible to have the LTI system for fixed P display a second-order 
response and it is important to note that this is not considered an 
oscillatory response for the purposes of this thesis. 

4.1 Linearized Analysis 

In Chapter 3, simulation results have indicated that the probability 
may tend to be neutrally stable in that it shows no tendency to return to 
its initial value following a perturbation. Greater understanding of the 
causes and consequences of this can be gained by examining the local 
behavior of the states and the probability about an equilibrium point. 
Examination of Equation (2.22) or (2.28) shows that the MMAC algorithm 
results m a closed loop system which has a set of equilibrium points; 
any point of the form w =0 , with any value of P is an equilibrium point 
since when w~0, there is no information m the system that would lead 
to a change in P.* The value of the state about which the present 
analysis is performed is w-0, P = l/2. The reasons for considering this 
point in preference to others stems largely from the high degree of 

symmetry inherent in the problem formulation given m Section 3.1. For 

\ 

example, the range of P for which the LTI system for fixed P is stable 


*Similar reasoning indicates that neutral stability is a general property 
of any adaptive controller. 
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is easily seen to be symmetric about P = 1/2 Ci.e,, about q=0). 

In order to lend clarity to the results, the continuous time equations 
are considered. Examination of the discrete time equations results in 
similar conclusions. The continuous time algorithm is rewritten as 

w(t) - A(q) w(t) (2.28) 

q (t) = ~(l-q 2 ) [r^r^ 1 jT 1 [ (l+q)^ + (l-q)r 2 ] 

= -Id-q 2 ) w' [£^ 0 . 2 ]w 
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where the linearization is done about the values [w^, q^] . Note that 

if w = 0 , then the linearized system reduces to 
”0 



Thus, the system is, at best, neutrally stable to first order (i.e., the 
linearized system has a zero eigenvalue) . This should not be particularly 
surprising; when the state is zero, there is no need or basis for changing 
the identification results. It should be further noted that if ACq^) 1S 
an unstable matrix, then the resulting linearized equations are unstable . 
Thus, m order to have asymtotic non-oscillatory behavior in a system sub- 
ject to perturbations, it is necessary that there exists a value of q Q such 
that A tqg) is a stable matrix. 

Equation (2.28) can also be rewritten exactly as: 


w(t) = ( 0) w(t) + Aw(t) q(t) (4.3a) 

q(t) = 7 w’ [Q 1 -Q 2 q - ^q 2 + O^lw . (4.3b) 

(This is the full Taylor series expansion about w=0, q=0) . A few ob- 
servations can now be made. First of all, for q near zero, the rate of 
change of the state is proportional to ] jw| j while the rate of change 
of the probability is proportional to | |w| | ^ - hence the neutral stabi- 
lity of the linearized system. Thus, the probability tends to change 
slowly for small values of w. Furthermore, again for q small, the higher 
order terms m q can be ignored so Equation (4.3b) can be approximated 


by 
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q(t) 




(4.4) 


Thus, the changes m q are such that q increases l-f | |-r- j | -> j l' r ^!'l _a nd 
decreases if | |_£ 2 j | < | |_r | | , which agrees well with the intuitive 
notions of how an adaptive controller should behave. Furthermore, when 
r^ = , there is no information as to which model provides the better 

match to the true system .and so no change in the identification results 
can be made. It should be noted that this may result in exponentially 
unstable behavior when Mq^-O) is an unstable matrix and the system is 
initialized with q(0) = q^ = 0 and, for example, r^(0) = £^(0) = _0. This 
is an important special case since in practice one often attempts to 
initialize an MMAC system with just such initial conditions. Equation 
(4.4) also illustrates the neutral stability of the probability. 

In summary, if MO) is an unstable matrix, the equilibrium w = 0, 
q=0 is unstable and small perturbations in the state cause a divergence 
from the equilibrium which will most likely result m the oscillatory 
behavior of the probability (resulting in either bounded or unstable 
state responses) observed in Figures 3.2 through 3.4. However, if A (0) 
is a stable matrix, then two modes of behavior are possible. For pertur- 
bations large enough that A(q) becomes an unstable matrix, at least 
temporary oscillations most likely occur as seen m Figure 3.6. For 
smaller perturbations, the state, w, will return to zero and the 
probability will simply move to a new value which is such that A(q) is 

stable. This occurs since, by Equation (4.4), the probability has no 

2 

tendency to return to zero unless | | | j - | | r^ [ [ 


changes sign. This 
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is precisely the behavior observed in Figure 3.7. 

The next section discusses the conditions such that the latter type 
of non-oscillatory behavior occurs for all initial conditions. This 
can be viewed as an extension of the results of this section for the 
case m which A(P) is stable for all P. Following an analysis of the 
stability of the oscillatory behavior, attention again is focused on 
the linearized system m Section 4.7 where a procedure is derived for 
determining how large the above perturbations can be and still have non- 
oscillatory behavior guaranteed. 


4.2 Universal Stability 

For any adaptive control algorithm one would at least like to be 
able to conclude that the overall system is asymptoticly stable about 
the point w=0 in spite of any uncertainties about the true system. As 
seen in the previous section, the MMAC method always results m a system 
with neutral stability in the probability. However, it is also shown 
that it is possible for the states to locally converge to zero despite 
the behavior of the probability. The present section discusses one 
case m which global asymptotic convergence of the states can be 
demonstrated . 

Consider the discrete time MMAC system given by Equation C2.22): 


w(k+l) = A(P) w(k) 

P(k) p^) 

P (k)p(r^) + (l-P(k)) p(r 7 ) 


( 2 . 22 ) 


P (k+1) 
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We now make the following assumption , which will be termed Universal 
Stability : assume j | a(P) | j £ e < 1 V P e [0,1] where j |a| | is the - 
hoffi of A. 'It should be noted that this is a- fairly restrictive 
assumption as it requires that the control law associated with each 
hypothesized model stabilize the true system. In general, this requires 
that the condition j | a(P) | | < 1 be explicitly tested for all values of 
P. However, as discussed in Section 4.3 in some cases such as when 
the true system and all of the hypothesized models are diagonal , it 
becomes sufficient to perform the test for extreme values of P only. 

When the assumption is valid, the following theorem is useful. 

Theorem : Under the assumption of universal stability, w(k) + 0 geo- 

metrically as k-* 30 . 

Proof: By assumption, | | | | £ £ < 1 *^£[0,1], Thn? p taking 

norms in Equation (.2,221 yields 

||w(jc+ll|| = j |a(p IwrCkl | | < ||a(p Lil ||wai|| 

— — 1 — — — i C4.5) 

£ e | |w Ckl | | 

Thus, j |w(k) | | £ e k | [wC.01 [ j , Since by assumption e< 1, the. conclusion 

follows . ^ 

It should be noted that by a suitable: redefinition of the. matrix 
2 ^ 0 ?} t this theorem can be extended to the N— model case where £ becomes a 
vector of probabilities such that | ]pJ | < 1 and each P^ is non-negative. 

The value of this theorem in applications is, of course, severely 
limited by the assumption on the stability of However, as 
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demonstrated by the analysis in the remainder of this chapter and the 
simulations of Section 3.2, a condition such as this xs necessary to 
guarantee geometric convergence, toy weakening of the hypothesis admits 
the possibility of at least transitory oscillatory responses. 

4.3 to Asymptotic Result 

As discussed in Chapter 1, the MMAC method was developed with the 
implicit assumption that one member of the set of hypothesized models 
exactly matches the true system. Under this assumption, Baram [15] has 
shown that the identification algorithm (i.e., the value of the proba- 
bility when all feedback gains are set to 0) converges to the matched 
model when the input is ergodic. However, he also shows that when de- 
terministic inputs are used and none of the models match the true system 
the convergence properties are indeterminate - the model to which the 
probability converges is a function of the input. Furthermore, his 
results require the ergodicity of the residual which clearly is not 
guaranteed when the probabilisticly weighted control is applied*. 

Thus , no general convergence result has been derived for the closed 
loop adaptive situation. 


*Ljung et al [24, 25] have considered closed loop identification for 
some specific model structures and identification methods. However, 
we have not attempted to apply their results to the MMAC method; 
this remains as a topic for the future. 
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However, the following theorem provides some insight into the 
expected behavior. It should be pointed out that no assumption as to 
the number of models is made. 


Theorem ; 

Cl) 

( 2 ) 

(3) 

(4) 

(5) 


Assume the following for a discrete tune MMAC System: 

There are N hypothesized models. 

Model #1 exactly matches the true system. 

All models (and therefore by (2) the true system) have 
diagonal A-ma trices with C=B=I (as in Chapter 3) . 

r^O) = 0 

P^O) > 0 


Then, x(k) 0 as 


Proof: Since _r^(0) = 0, and model 1 matches the true system, r^(k) =0 ¥k. 

Thus, P^(k) is non-decreasing and either: 

Case 1) P^(k) ->1 as k-*- 03 or 

Case 2) ?e>0 3 1-P (k) > e Vk. 

(equivalent, 3i3P 00 > 0 .) 

Consider Case 1 : 

P^(k) +1^3 K 3 k > K — >A(P (k) ) is a stable matrix since A_(l) is 
stable by hypothesis . This m turn implies the convergence of | j x [ ] . 

Consider Case 2 : 

l-P(k) > £ Vk together with P(k) non-decreasing implies 9- at 
least one other model such that its residual approaches zero since 
otherwise P^(k) would increase to 1. Assume for the moment that there 
is only one such other model, denoted by model i. (Thar is, assume 
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that only r., CO and r. approach- zero,) Due to the use of full state 

— T_ x ■’ ■ ■ ■ * • • 

observations in the canonical problem , the equation for r^ CO can be 
rewritten from Equation C2.22) as 


r (k+1) » A (IOH) r (10 + (A-AjxCk) 

—1 —i -i i — 


Since r^ (k) -+ 0, ic (3c) must approach, an element of N(A-*AO,* The dynamic 
equation for x(k) is now rewritten as 

N 

x (k+1) = B (k)x(k) + E P. (k)G. (I-H.)r. (k) 

-0 - j=1 3 "3 “3 -3 

where 


B-(k) = [P.(k) (A-G) + (1-P. (k) ) (A-A ) + (1-P <k) ) (A -G. ) 

—0 1 1 1 1 —1 —i 

N 

- E p (k)G ] . 

3 =2 3 "I 


If r.(0, j^l or i is bounded, then, since P . (k) -+ 0, P G (I-K.)r.(t) 0. 

“3 J 3 3 3 3 3 

If r . ( O , jrl or V 03 ' then, since p, +0 as e~ r 3 r 3, p. r + 0. 

~3 3 3 “3 

N 

This implies that I P (k)<3. (I-H )r (k) -*■ 0_* Thus, consider the un- 

j = l 3 3 33 

undriven system 


x(k-KL) = S Q Ck)xCk) 

which, since P^C 00 ), P (“) exist, can be rewritten as 


*N C* ) represents the null space of a matrix 
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xOc+1) = [P x (<»} (A-G^) + Cl-P-^ 00 )) (A.-G.)]x(k) 

+ [(P 1 (k) -P^)) (A— G^) - CP 1 (k)-P 1 (<»)) (A^G^JxOt) 

+ (1-P^ (k) ) (A-A^)x(k) 

(A-G^) and (A^-G^) are stable by design. Then since they are also 
diagonal, the convex combination of stable matrices is stable. Thus, 
there exists a b<l 3 

i | [P <~) (A-G^ + (1-P (*) ) (A^-G^ ) ]x(k) | | £ b| |x(k)| j 

Then , for some £ 3 b + E < 1 and any given 6 , we can find a K 3- 

j | [(P <k)-P (»>) (A-£ )-<P (k)-P (»)) (A i -G i ) ]x(k) | [ £ £ | [x(k) | | 

(since P^(k) ■+ P^ 00 )) and 

| j (1^00) (A-Ajxtk) | j £ 6 
(since x(k) + N(A-A^)). Combining these yields 

| |x(k+l) | | £ (b+£) |-|x(k) j | + <5 ¥k > K 

or, using the variation of constant formula, 

k-1 

| [x(k) I I £ (b+£) k “ K [ |x(K) ||+ 2 (b+£) k ~ 1_:, 6 ¥k > K 

j=K 

or | |x(k) j [ £ (b+£) k K j |x(K) | | + ^ $ — 


Taking the limit as k->“ 
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Lxn | jxOc) | | £ 0 + 1 

since K was selected such that tb+e> <1. But S' was arbitrary. Thus 
Lim | |.x (k) || =0. The case m which more than two models have nor.- 

JC+OD 

vanishing probabilities follows similarly. 

A 

Some comments on this theorem are important. 

1. Assumption (.5) merely guarantees that the probability of model 1 
can change since 

p CO) - o => p ( k) =o vk . 

i 

2. Assumption C4) has been made for convenience. It is used to 
guarantee that ( - ) is non-decreasing, thus preventing oscillations in 
the probability due to | l^t*) j | - J (r^t*) | | changing sign. 

3. Assumption C3) is, of course, very restrictive. It is necessary 
to enable one to conclude stability when a convex combination of stabi- 
lizing controls is applied. Given a non-diagonal matrix A and controllers 

and with 

IlCA-Gj | < 1 , i = l, 2 

it can be shown by counterexample that 

A-bG -(l-b)G be [0,1] C4.6} 

- J. 4 

is not necessarily a stable matrix. It may be possible to overcome this 
problem by using the fact that Ck) 0 as h + 40 which restricts the 
kinds of state interactions which can occur in Expression (4,6}. It 
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should be pointed out that assumption C3) can he relaxed to any situation 
in which the convex combination of stable matrices is stable. For example, 
if A and all of the A are in Jordan form with identical Jordan structure, 

— —i 

then our result holds. This is equivalent to assuming that the shapes 
of the modes of the system are known but that the eigenvalues are not. 

This is less restrictive than the diagonal case. 

4. Note that no claim is made about the convergence of the 
probability even when one model matches the true system. This is, of 
course, another example of the need for persistently exciting inputs 
in order to guarantee the convergence of identification algorithms [36] . 

4.4 An Analysis of the Probability Equation 

The simulations of Section 3.2 have indicated that the behavior of 
the probability often resembles the output of a switch, alternately 
taking on values near zero and near one but seldom in between. By now, 
it should be clear that this property is largely determined by the 
equation for the probability. Equation C2.22b). Therefore, the present 
section contains a detailed examination of the characteristics of this 
equation. The principal conclusion is that the equation for the 
probability can be rewritten as a scalar, static nonlinearity and 
a summation (.the log likelihood ratio) . This decomposition aids the 
analysis since attention can be focused on each part separately . Thus, 
this section examines the characteristics of the nonlinearity and shows 
that the switch-like behavior of the probability is largely due to this 
nonlinearity. Section 4.5 continues the analysis by examining m de- 
tail how the true system and set of models affect the log likelihood 
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ratio and therefore, by the analysis of this section, the probability 
and ultimately the closed loop behavior. 

Equation (2.22b) is repeated for convenience: 


P (k+1) = 


where p(r^) = Be 


P(k)p(r 1 ) 

P(k)p(r 1 ) + (l-P(k) )p(r 2 ) 

This can be rewritten as 


(2.22b) 


P (k+1) 


P(k)y(k+1) 

P(k)Y(k+l) + (1-P (k) ) 


(4.7) 


where Y(k) is the instantaneous likelihood ratio 


^[r-0- 1 


Y(k) = 


p (r x (k) ) B 1 e 
P(r 2 (k) ) = “ 


2 till ^l 3 


3 e 
P 2 


-— [r'0 1 r ] 
2 — 2—2 —2 


= B e 


- — fr' 0 -1 r - r ' 0 _ ^r 1 

2 — 1— i -i —2~2 — 2 J 


(4.8) 


It is now possible to rewrite Equation (4.7) such that the probability 
does not appear recursively: 


P(k) = 


p (o) n y (k) 

1=1 


p (o) n Y(k) + (i-p(o)) 

1=1 


(4.9) 


Finally , let 


k 

E 

i=l 


Fr 1 0 ^"r 
L -l-l -l 


r ' 8 ^r ] 
— 2—2 —2 


a(k) 


(4.10) 
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Thus , Equation (4.9) becomes 


„ - |ct(k) 

P(0)6e 

P(k) - (4.11) 

P (0) Be + (1-P(0)) 

Note that ct(k) is the Log Likelihood Ratio. Thus, Equation (4.11) 

provides the connection between the Log Likelihood Ratio and the 

probability. Figure 4.1a is a plot of P Ck) versus a(k) for a few 

values of P(0) and {3 = 1. Figure 4.1b is a detail of the same curve 

for P(0) = 1/2. Define a to be the value of a(*) for which P(-) = 1/2. 
s 

It ■ is then clear that 


a(-) 

a(-) 

Thus , a will 
s 

can be solved 


a = 
s 


» a => P ( • ) « 0 
s 


« a => P { • ) « 1 * 

s 


be called the switch point for a ( • ) 
to give 


-2 In 


E o B 


(4.12) 


Equation (4.11) 


(4.13) 


Now, consider Equation (4.11) evaluated m the vicinity of a . For 
example , if 


a( •) = a + 2e 
s 

where e is any positive or negative number we see that the resulting 
probability is 

-e 
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a(k) 

) versus a ( • ) 


of = 1/2 
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whxch is totally mdpendent of any parameters of the MMAC method. Finally, 
simple calculations reveal that the curve of P ( • ) .versus- a(‘) is anti- 
symmetric about the point a . Thus we see that the switch point a is 

s s 

determined solely by the aprion information P(0) and (3 and that no 
parameters of either the true system or the models (except for the 
assumed noise covariances contained in $) affect the switching character- 
istics of the system other than by determining ct(k) . Various approxi- 
mations to P{*) can thus be calculated from a(*). One example , shown 
is the dotted line in Figure 4.1b, is 


P (k> = 


1 

0 


a( • ) < a s - 5 
a( •) > a +5 


(4.14) 


s 

\. 5-(a(-) - a )/10 otherwise. 

' s 

Examination of Equation (4.11) or Figure 4.1 reveals the reason for 
the zero-one type behavior so often noted in the simulation results m 
Chapter 3. Only for values of ct ( • ) such that 


|a(-) - a g i < 5 

will an intermediate value of P occur. As seen m Equation (4.10) 

(see also Section 4.1), Ct(0 is proportional to the square of the norm 
of r ( • ) . Thus , for r ( • ) large , a ( * ) tends to change at a rate twice 

—i —i 

as large as that for | |r^(>) | j and may never fall m the range. 

One potential advantage of using Equation (4.11) m a MMAC con- 
troller is that a wealth of information has been accumulated about 
the behavior of the log likelihood ratio [34, 35]. Thus, use of this 
approach allows full use of this information while still permitting a 
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sxmple control to be calculated. 

The approximation (.4.14) is. also useful for characterizing the . 

length of the half-periods of the oscillations observed in Section 3.2. 

By the definition of a , it is clear that [jumps in the probability (i.e., 

s 

from 0 to 1 or from 1 to 0) correspond to transitions of a ( • ) through 
a^. It is thus possible to calculate a bound on the half -period of an 
oscillation. Define T^ to be the time of a transition m the probability 
(assume from P = 0 to P = 1) and (T + T^) to be the next transition. (see 

Figure 4.2). Thus, by the above analysis and the definition of ct s 

a(T n ) = a = a(T. + Tj . (4.15) 

Is 12 


By the definition of a(-) 


£ (i) 5 X £.i ti) ~,r^(i)^ X r (i)] 

i=l 


T 2 +T l 


= £ [£j_(i) i^d) - (DS^r. 2 (i) ] 

i=l ~ 


■ I ^ A -1 > 


T +T 

21 , , 

° r £ (r ' (i) e. V (i) - r ' (i) 9~ r (±) ) = 

i=T +1 _1 _1 _1 ” 2 “ 2 “ 2 


(.4.16) 


C4.171 


Equation (4.17) thus provides a condition which. and must satisfy. 
This equation is ejqplored further in the next section where an approxi- 
mation for r. (*•) is employed. 

— i 

This section has presented a detailed analysis of the probability 
equation (2.22). It has been shown that -the equation can be broken 
down into a static nonlinearity and a summation (.the log likelihood 
ratio) . Attention has then been focused on the static nonlinearity. 
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Fig. 4.2 Switching Times 
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The next section continues the analysis by considering the hehavior of 
a(.) for one class of systems. 

4.5 An Analysis of the Oscillatory Behavior 

Many of the simulations of Section 3.2 display an oscillatory be- 
havior involving all states and the probability. The preceeding section 
contains an analysis of the equation for the probability of a model, one 
of the conclusions of which is an observation that Equation ( 2 . 22 b) can be 
divided into two parts; one containing a static nonlinearity C4.10) and 
a relatively simple summation 


a(k) = Z ||r (k)H 2 - !|r (k)|| 2 

i=i 1 q/ 2 e/ 

—1 —2 

(4 *15). 

k -1 -1 
= Z - £ 2 (k)S 2 £ 2 DO 

1=1 


Furthermore, the approximations 


a(-) » a ?=> PtO » 0 
s 


C4,m 


at*) « a => PC*) *• 1 
s 

have been found to. be very good for [a(-) - a | greater than about 5 

S’ 

where a^ can be determined' solely in terms of PCO) and 6, Note that 
for the present study we can take a = 0 since we can assume that the 
system started in the remote past in which case P (0) can be chosen 
arbitrarily . 
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The present section continues the analysis of the. oscillatory hehavior 
of the probability by considering the behavior of the variable a(k) for 
the canonical problem of Section 3.1. The approach, taken is to isolate 
the basic modes of the state of the canonical problem by partitioning w 
m a new fashion. This partitioning emphasizes the two basic modes: 
inherent in the structure, namely those associated with true state #1 
and those associated with true state #2., Because of the diagonal nature 
of the canonical system, this partitioning allows aC*) to be bounded by 
simple exponentials which are then analyzed. 

In this section extensive use has been made of the observation 
from Chapter 3 that during periods of oscillatory behavior the prob- 
ability tends to be virtually constant for long periods of time and then 
abruptly changes. This square-wave like behavior, clearly seen m 
for example Figure 3.2, is a key element of the approximations employed 
in this section. 

Recall from Chapter 2 that 

wCk+1) ~ A{P)wCk) C2.221' 

where w'CO = [x 1 (.* 1 (., ) r | (_• ) ] . As seen in Section 3,2, during 

periods of oscillatory behavior, each component of wC* } is alternately 
stable and then unstable. It thus is natural to regroup the states 
such that components which are simultaneously stable are grouped to-* 
gether. To this end, define y^ C * 1 to be those states of w which cor* 
respond to the first component of the true state and y^ C« ) those which 
correspond to the second. That is, 
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*1 c-r 


l 

/•“N 

« 

CN 

£ 



; y 2 XO = 

r^ 2) (•) 




rfh.) 


where r^^ C*) is the x^- component of the residual of the KF for the 
j— model and x^ ( • ) is the 1 — — component of the true state x. Using 
this decomposition Equation (2.22) can be rewrxtten as 


'y^k+1) - 


^(P) 

0 


i 

N 

H 
< — - 



£ 2 (k+1) 


J3 

V p) - 


U 2 {k) J 


(4.18) 


where A^(P) and (P) contaxn the appropriate elements of A(P) . 
block diagonal nature of Equation (4.18) is a direct consequence 
the assumption that the true system and each model is diagonal. 

and are defined to be the appropriate partitioned versions 
j) ^ , Equation (4.15) becomes 


The 

of 

If 

of 6_ ^ and 


L k i l 

2 y^(0) nA^(P(D))^ 2 H a (P(3))y (0) - 

=1 3=1 3=1 

^(o) ns i (p(3))$ 1 n a 1 (p( ] ))y 1 (o) . 

1 3=1 3=1 -* 


(4.19) 


We will consider the initial conditions shown in Figure 4.3(a) 
in which a(0) < 0 (but a(-l) > 0) , | jy^O) | | 2 > | |y 2 (0) I 1 2 and P = 1. 
Note that this corresponds to the case m which a change in probability 
(from P = 0 to P = l) has 3 ust occurred. The equation 04 .19.) and the- 
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simulations of Section 3.2 indicate that if oscillatory behavior occurs, 

the traj ectories of | | y C * ) J | ^ , a ( • ) and PC*) must be as sketched in 

Figure 4.3, since [ [y (.0 [ [ is decreasing by design for PC*) ” 1- 

It will be assumed for the present analysis that P can be closely 

approximated by either 0 or 1. As seen m the previous section, this 

corresponds to assuming that aC*) - a I is greater than about 5. From 

s 

the simulations of Chapter 3, this is a reasonable approximation during 
periods of oscillatory behavior such as are being analyzed here . * Further- 
more, it will be assumed that both K.F. ’s have been correctly initialized 
(equivalently, have been running since the remote past) so that the 
matched K.F. residuals are zero. 

The remainder of this section is concerned with characterizing the 

peak excursions of J |y |( (see Figure 4.3) and also the half-period T^. 

This is done by approximating the behavior of a(*) after first approxi- 

1 k ~k 

mating the behaviors of y' (0)A. (b.A (l)y (0) (see Equation 4.19). It 

— l —i —i—i —i 

is clear that each of these terms must behave as exponentials for P = l; 
that is 

||y i (k)|| 2 < aj ||y i (0)|| 2 (4.20) 

2 

where a^ = j | a^( 1) | j . Thus 

« 

a(k) « E [a^ ||y 2 (0)|| 2 cJ 2 -a^ |iy 1 (0)|| 2 a i ] (4.21) 

i=0 

* 

When the probability is forced to either zero or one Cby, for example, 
use of a maximum likelihood type control rather than the probabilistically 
weighted one) then this is not an approximation. 
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where = the maximum eigenvalue of 9 ^,* A few comments are necessary 
on this approximation. First of all, note that for P = 1 it can be shown 
that 


a 


a 


1 

2 


< 1 
> 1 


by the constructions of A ^ ( • ) . Furthermore, as noted earlier, 


llx 2 (o)!l 2 « IIy^o )!! 2 • 


14.22) 


Thus, during the initial part of the interval, the term m y^ ( • ) dominates 
while later the term m y (•) dominates. We are interested m determining 
the half -period, i.e., when the probability switches from 1 to 0. As seen 
in Section 4.4, this occurs when Ct ( * ) =0. Thus, an approximation to 
this condition is to find such that . 

T 1 . T 1 . 

S S lly?(0> I I 2 * S ! 1 7-i (0) I | 2 . (4.23) 

i =0 A * 1=0 1 

It should be pointed out that the term in y^ leads to an underestimation 
of T 1 while the term in y. leads to an overestimation. Thus, Equation 
(4.23) does not yield either an upper or lower bound. However, it is 
felt that this is still of value since the analysis yields insight 
into the type of behavior observed m simulations. Furthermore, al- 
though sometimes not of use, the of Equation (4.23) does lie between 
the upper and lower bound. 


*Ry the symmetry of the canonical problem, note that = (T„. 


Evaluating the summations m Equation (.4. 23) yields 


- il 


fa' 1 !* 11 „ ll 

1 

15 liy 2 (o)II 2 

fM 

a. - 1 

a. - 1 

1 


M 

1 


where by the assumption that on the interval k=0 to T^, P « 1, 

= ! 1 ^( 1 ) | | < 1 

* 2 = I|a 2 u)M > 1 


and use has been made of the fact that CT - L " 0 2 * Equation (4.24) could 
be solved numerically to arrive at an approxxmati on of although no 
closed-form solution is possible. 

Further insight into the behavior can be obtained by invoking 
another assumption. Since 


and 


< 1 
> 1 


we assume 


T +1 
a l 


« a. 


V 1 

a 

2 


» a. 


(4.25) 


This is clearly reasonable since T^ depends on a^ and a^ m such a 
way that a^ closer to 1 yields larger. Thus, although (4.25) is 
an assumption, it is very reasonable. Using this. Equation (4.24) 


becomes 
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ll (0) 


(1 " a l ) 


Y 0 CO) | 


T +1 

_ a 
2 2 


a 2" 1 


or, using Equation (4.20) 


II 2 * I 1^(0) | i 2 (&2 


a 2 (l - a x ) 


Defining 


R = 


(a 2 " 15 
a 2 _(1 " a i } 


(4.26) 


(4.27) 


(4.28) 


yields 

IIy 2 (t 1 )!I 2 * R l lx 1 (°> 1 1 2 • 


(4,29) 


Thus, if R > 1, the value of j |y^(«) | [ when P switches increases with 
each cycle and ] |y {•} | | -*• 00 (i.e., the closed loop system is unstable) . 
Similarly, if R < 1, then peaks of | |y^(») 1 | decrease with, each cycle. 

' An interesting interpretation of these results can now be given by 
noting that 


a 2 ^ 1 — — ^ R > 1 (4.30a) 

a 1 a 2 <1 ^=>R< 1 (4,30b) 

a i a 2 = 1 ^ R = 1 (4.30c) 

and that a^ is the slowest decaying mode of the stable matrix A^ (1) and 

a 2 is the fastest growing (most unstable) mode of A 2 (l). Thus the 
closed loop system will be stable if the slowest decaying stable mode 
of the matching subsystem (A^ (1) ) dominates the unstable modes of the 
mismatched system m the sense that the product of the appropriate norms 


is less than 1. 
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Equation C4.26) can also be rewritten as 

i-2 


T x III! CO) I 


a 2 1 


|y 2 (°) 


R 


£4,311 


Thus, the following approximation to can be given 


T x — In 


III! C0)| 

2 

1 Iy 2 co) | 

2 


/ In a. 


£4,32) 


An interesting phenomena can be seen by considering multiple periods. 
Thus, by analogy with Equation (4.32) 


T 2 «« In 


R ! |y 2 (T i J I 

II I 

But it can be shown that 


/ In a 2 


l|y 2 (yll 2 S 1 1 y x CO) 1 1 2 

II^CTi)!! 2 a ^i 1 1 y. 2 Co ) 1 1 2 

This results m the conclusion that if the initial period is sufficiently 
long, then subsequent periods will be of increasing length if a^a 2 < 1. 

The above analysis results m the conclusions that a ]_ a 2 t * ie 
key quantity m determining the stability of the closed loop system. 

If R is less than one, then the peaks of the state trajectories decrease 
and the period increases. However, if R is greater than one, instability 
results. These conclusions agree well with the simulation results of 
Chapter 3. The only point of discrepancy is in the simulation of Figure 
3.2. In that case, the peaks appear to remain constant even though 
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calculations reveal that stability results. Close examination of the results 
of the simulation have revealed that in this case the lower limit on 
the probability (see Chapter 3) is consistently achieved. As shown in 
Chapter 5 , this somewhat modifies the behavior resulting in uniform peaks 
even when the system is stable. Further discussion of this important 
case is deferred to Chapter 5 . 

In summary, the stability of the MMAC algorithm is determined by the 
relation between the growth rate of the most unstable mode of the mis- 
matched subsystem and the rate of decay of the slowest stable mode of 
the matched subsystem. This leads to the qualitative conclusion that 
stability results when the stable modes are "faster" than the unstable 
modes. This concept is further emphasized m the following section 
where an alternate view of the problem leads to a similar result. 

4.6 Quasi-Lyapunov Analysis 

An alternate view of the stability results given m the previous 
section can be obtained by employing an analysis which closely resembles 
the methods of Lyapunov. The results obtained compliment those of 
Section 4.5 and given increased insight into the behavior of the MMAC 
method. 

The approach of this section is first to eiqplore the use of the 
"normal" quadratic Lyapunov Function. This is used to introduce the 
ideas and also to demonstrate the drawbacks of this standard approach. 

This leads to the introduction and motivation of a new Lyapunov-like 
function which is then used to derived a stability result which. 
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emphasizes the role of the individual control gains. This function is 
motivated largely by the simulation results in which the controller 
is seen to alternately attempt to control each mode of the system. Thus, 
each state tends to be piece-wise exponential Csee Chapter 3) . This 
results m a phase-plane plot, such as seen in Figure 3.5 which resembles 
a hyperbola. It is this observation which led to the function investigated 
m the major part of this section. 

First, consider the function 

v(k) = w' (k) w(k) C4.33) 

where w(k) is governed by Equation C2.22) 

* 

w(k+l) - A(P) wCk) . C2.22). 

Thus, vtk+1) can be written as ' 

v(k+l) = w 1 tk) A 1 (P Ck) )_A (P Ck) )w(k) . C4.34) 

In order to prove stability, it is required to show that 

v(k) > vtk+1) . (4.35) 

Unfortunately, Equation L4.34) depends- strongly on P and, thus stability 
is difficult to show. The one case in which it is possible to determine 
stability is when | |aCP) | | <1 VP, However, this case has been well 
discussed in Section 4.2 and so will not be repeated. 

An improved result is obtained if attention is first restricted to 
the diagonal case. Thus, assume that the canonical system of Section 
3.1 is under consideration. Define the function 
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v(k) = x^(k)x (k) 


(4.36) 


or 


v(k) = — w' (k)Iw(k) 


(4.37) 


where 


I = 


"0 1 

1 0 


A definition is useful in the sequel: The system (2.22) is hyperbolically 
stable if v(k) 0. The choice of v(k) is heavily motivated by the 
simulation results in Chapter 3. As discussed in Section 3.2, the phase- 
plane plot of vs. x 2 resembles a hyperbola (truncated for large 
excursions). Figure 3.5 graphically depicts this phenomenon. A second 
motivation stems from the analysis in the’ previous sections of this 
chapter in which each state tends to be governed by an equation such as 


x. (k+1) 
i 


= a k x. (k) 
i 


(4.38) 


for periods of time. This, of course, is due to the fact that the 
probability tends to be piece-wise constant due to the nonlinearities 
analyzed m Section 4.4. 

It should be pointed out that the function in Equation (4.36) is 
not a valid Lyapunov function. First of all, v(k) can be either positive 
or negative depending on the initial conditions. Secondly, v(k) = 0 
does not imply that x^ is even bounded 1 Neither of these invalidate 
the analysis, however. A slight redefinition of v(-) would overcome 
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the first. This is not e:xplicitly done here because it is not necessary 
for the problem at hand and would complicate the analysis. The bounding 
of the state- has al-ready been indicated in previous sections and thus 
is provided by alternate means . Furthermore , even if a bound were not 
available , the analysis would yield insight into how the system behaves 
when each mode is alternately controlled. 

Some useful notation for the diagonal system is now needed. Recall 
from Section 3.1 that 


Define 


and 



"a 0~ 


a 0 

A = 


ii 



0 a_ 


0 a. 


’g O' 


0" 

V 


; —2 " 



[0 gJ 


.0 g_ 


■h o' 


“fi o' 

H- 


t h = 


— X 

_P 

z 

-0 h. 


Temporarily, assume that g = 0.* Then Equation (4.37) becomes 


v(k+l) = w' (k)A’ (P)J«(P)w(k) (4.39) 

which, when written out term by term becomes 

2 2 2 

v(k+l) = [a(a-g) + Pg -P g ] v(k) (4.40) 

or v(k+l) = A(P) v(k) . Note that A(P) , which is plotted in Figure 4.4, 


This corresponds to Case 1 in Chapter 3. 
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is a convex quadratic function of P. The following two theorems then 
follow : 

Theorem ; Given the canonical MMAC system of Section 3 . 1 with g = 0 , 
the overall system is hyperbolically stable (i.e., x^x^ 0) if the 
following two conditions are met. 


(1) 

|a(a-g) | < 1' 


(4.41) 

and 

(2) 

Ja(a-g) + - g 2 | < 1 

* 

(4.42) 

This plus the bounding results 

of Section 4.5 yields 

stability. 

Theorem; 

Given the canonical MMAC system of Section 

3 . 1 with g =* 0 , 

the overall system is unstable 

if 


(1) 

a ( a-g) > 1 


(4.43) 

or if 




(2) 

/ » 12 
a(a-g) + -j g < -1 . 


(4.44) 


A few comments can now be made. 

1. Conditions (4.41) and (4.42) essentially require that v(k) 
be decreasing for all p . Note that this can be -true independent of 
whether A(P) is stable. In fact, for the Case 1 example of Chapter 3, 
A (P) is instable VP but v(k) is still decreasing. 

2, Theorems such, as the above can easily be derived for the 
relaxed case in which the two elements of the true system matrix are 
not equal. The curve for A(P) is -then not symmetric about P = 1/2 but 
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similar analysis can still be done. The results do not significantly 
increase our understanding and so are not presented. 

(.3) Note that Condition C4.431 is consistent with the results sum- 
marized by Equation (4.30) of the previous section. The added condition 
of C4.42) results from relaxing the assumption m Section 4.5 that P 
is always "near" zero or one. Further insight results from plotting the 
Conditions (4.41) and C4.42), The shaded area in Figure 4.5 is the 
set of a and g for which Condition C4.41) is satisfied and the crosshatched 
area is where Condition (4.42) is satisfied. The points corresponding 
to Cases la, b and c from Section 3.2 are plotted. 

A few comments are in order. First, as evidenced by the simulation, 
results , it may be sufficient for hyperbolic stability that Condition 
(4.41) be satisfied. In that case, if p oscillates between zero and one 
(which the simulations indicate is likely) , stability will result. 

Further, if the probability is forced to be either zero or one by using 
the control for the most likely model rather than the probabilistically 
weighted control, then Condition (4.411 is sufficient to guarantee 
stability. 

Another interesting observation is that for this special case, the 
smallest control gain which yields vC-) decreasing for p = 0 or 1 (i.e. 
the lower edge of the shaded area in Figure 4,51 exactly corresponds to 
the value of the control gains for the LQ problem, with a state cost of 
zero I 31 , 32 ] . This implies that for this special case use of the LQ 
methodology guarantees that Condition (4,41) will be satisfied. This 


is believed to be due to the robustness properties of the LQ design 



106- 



Fig. 4.5 Regions of Decreasing v(-) 
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technique discussed in [5 , 33] , The implications to more general MMAC 
structures remains an open question. However, note that there are many 
control gains ( not generated by an LQ methodology) which yield unstable 
behavior. For example, the matched system is suable for all control 
gains 

g e (1,2) 

However, from Figure 4.5, if the chosen control gain is 
g e (1, 1.5) 


then the resulting MMAC system is unstable. 

A second point to note from Figure 4,5 is that for a >2, no linear 
controller will result in vtk) decreasing for-P = 1/2, Under these 
conditions , the only stable behavior possible .must be of the oscillatory 
type. 

The removal of the assumption that g = 0 results in the analysis he^ 
coming intractable. For example, when evaluated at 5 = 1, the equation 
corresponding to (4.40) becomes 


v(h+l) = (a-g) (a-g) v(Ic) + 
+ (a-g) g(l-h)r^ 1) x 2 

+ (a-g) gCr^r^x^ 


(4.45) 


+ g 2 ( 1 - h ) 2 


r 


( 2 ) ( 1 ) 

1 r l 


If it is assumed that the matched KF state has been correctly initialized 
(i.e., r^"^ = 0) this becomes 
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v(k+l) = (a-g) (a-g)v(k) 

+ (a- g)g(l - fijr^x^ 


(4.46) 


where the mismatched component of the KF now enters to modify v(k+l). 
Equation (4.46) has evaded further analysis except in a few obvious 
special cases (such as 6=1) which will not be discussed. 

This section has attempted to accomplish two goals . First of all , 
one interesting special case has been examined to the point that both 
stability and instability results have been presented. These results 
complement those of the previous section and provide an alternate view 
of the mechanisms underlying the stability of the MMAC method. A pos- 
sibly more valuable contribution has been to point out the need for 
careful design of each 'individual controller and the set of available 
models. 


4.7 Domain of Attraction 

In Section 4.1, it is demonstrated that if A (p) is a stable matrix 
for P = 1/2 , then the system linearized about P = 1/2 is neutrally stable 
in that the probability has no tendency to return to 1/2. This is seen 
to be due to the fact that if w = 0 , then there is neither the need nor 
the basis for changing the probability. Equivalently, the equilibrium 
set is {w = 0 with any p}. This is sketched in Figure 4.6. The stability 
of each equilibrium point has already been examined using a linearization 
approach. One important question which remains, is to determine the set 
of initial conditions such that the system will return to an equilibrium 
point without first oscillating. (The analysis of Section 4.4 through 4.6 
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have indicated that stable oscillations may result in w approaching the 
equilibrium set and it should be noted that this is not the case under 
consideration here) . 

The importance of this analysis is best seen by considering the 
simulations of Case 2 in Chapter 3. There it is seen that if A(P) is a 
stable matrix for P=l/2 but not for p = 0 or 1 then small initial con- 
dition perturbations result in small changes in w ( - ) and P(-) whereas 
large perturbations result in P ( * ) going to zero or one and then os- 
cillating at least for a while. This section continues the analysis 
of Section 4.1 and results in a procedure for determining a bound on 
[ Iw(O) | | such that the probability does not oscillate. 

Consider the equations from Section 3.6.1 


w(k+l) = A(P) w(k) 

P(k+1) = P(k) p(£ 1 ) 

P(k)p(r 1 ) + (l-P(k) )p (r 2 ) . 


( 2 . 22 ) 


We will assume throughout most of this section that P(0) = 1/2 (i.e., 

a =0)* and that A{*) is such that there exists an £ < 1/2 such that 
s 


| | A(P) | j < 1 V P £ a/2 - £, 1/2 + £) 

1 I^CP) I I >_ 1 otherwise. 


C4.47) 


It can be shown that this is the only case in which exponential behavior 
has not been investigated. In Section 4.2, the case of j |a(P) | j <1 VP 
has been investigated, and if | |a(P) | | > 1 VP then Section 4.1 has- 


*As m Section 4.4 we assume 3=1- 
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shown that oscillatory behavior (either stable or unstable) must result. 

Assumption (4.47) represents the examination of a previously unstudied 

case. As seen in previous sections, A(P{-)) is a function of ci(0 which 

in turn is a function of the data (state) generated by A(P ( - ) ) . Thus, 

there is a close connection between bounds on ct(*) and bounds on 

| |a(P(*)) | |. The basic approach of this section is to bound «(*) in 

such a way that | |a(P( 4 )) ] | remains a stable matrix. This bound on 

a(‘) is then translated to a bound on | jw^ | |- That is, we bound the 

size of the initial conditions so that w(k) is small enough so that 

-0 — 

in turn a(*) remains within its bound. This is accomplished as follows. 
Prom Equation (4.48) it clear that 


' p e (~ - e, J + e) => | |£{p) 1 1 < l (4.48) 

and that P £ - £., , j + => | |a(P) | | < a < 1 for some ^ < £. 

(Refer to Figure 4.7.) From Section 4.4, each value of P(*) maps to a 

value of a(-) as shown in. Figure 4.7 and thus attention can be focused 

on a(*) . Define 6 to be the value of a(*) such that P = 1/2 - £ and <$^ 

as the value of a(*) for P = 1/2 - £^. 

It is shown in Appendix A that 
ii i i2 


a(k) 


*0 I 


- ( 1-^)0 


Vk 


(4.49) 


where a is Max | j_A(P(k)) | j* and a is related to the residual weighting 
1 a(-) 

matrices of 0, and EL . Thus, one may choose an a{») =6. for some 6 , 
—1 -2 11 


*The maximum is achieved for the maximum value of a(*) for the canonical 
problem. 
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compute the value of a < 1, and finally use the equality from Equation 
(4.49) to find a | j | ] 2 such that a(k) is bounded as necessary. By 

2 

changing 8^ and repeating the procedure, a different value of j [ w^ [ | 
is obtained which may be either greater or less than the first. The 
detailed development of the algorithm is gxven in Appendix A. The re- 
sulting procedure may be formalized as follows : 

Step 1 ; Find e such that Condition (4.48) is satisfied. This in general 
will require an iterative procedure for determining [ ]a(P) | | for various 
values of p. 


Step 2 : Calculate 6 from Equation (A.l) . 

Step 3 : Choose an 0 < < 8. 

Step 4 : Calculate corresponding to from Equation (4.11) . Finally, 

compute a = | | A(P^) | | . Note that by the selection of a in Step 3, a 
is less than 1. 

2 

Step 5 ; Calculate | jw| | from 


m2 aa(l-a) 
„|| = 


( 4 . 50 ) 


Repeat steps 3 through 5 for different values of to maximize j |wj | . 

Note that a maximum exists since all functions used above are continuous 
and | jw| j is minimized for a — 0 and for a =-6. The best way to understand 


. 2 

Due to the nature of the canonical problem, | ]a(P) | [ is maximized for 

extreme values of P. In general, greater care is needed in order to 

2 

assure | |a(P) [ | a for all p in the region of interest. 
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the procedure is to verify that [ |w| | in Equation (4.50) results in 
exponential behavior. To do this it is sufficient to show that 
I ! w 0 I I = 1 l w l I =>|cc(k) | < <5 Vk . But by construction 

|a(k) | £ a Vk 

and so a < 6 by Step 3 finishes the result for = 1/2 . 

As an example of the use of this algorithm, Table 4.1 demonstrates 

its application to the conditions of the Case 2 simulations from Chapter 

2 

3. The table contains the evaluation of [ [w_ | | for various values of 

It is thus possible to conclude that if P(0) = 1/2 and | |w(0) | | £ .236 

then exponential behavior will result. Simulations confirm this obser- 

2 

vation and indicate that | jw(0) | j .8 results in P just reaching the 
boundary of the stability range of A(p) . 

It is possible to extend the procedure to P^ S (1/2 - £ 1/2 + £) as 
follows. Given p , compute a Q from Equation (4.11). This effectively 
reduces the amount a(*) can change and still have A.(P(»)) stable. 

Thus replace 6 in Steps 2 through 5 by 5- |a Q |. It thus becomes pos- 
sible to compute a complete region in P^ x ] [w^ | | space such that a 
nonoscillatory response results. Figure 4.8 depicts such a region. It 
should be pointed out that the resulting set is not exhaustive in that 
points outside the set may result in nonoscillatory behavior. This is 
due largely to the conservative nature of Inequality (4.49) which uses 
a worst case estimate of | jA(P(k)) | | for all k. Thus, a(k) tends to 
change slower than predicted by Inequality (4.49). The procedure does, 
however, yield a useful lower bound to the full set of mitral condi trons 


which yield nonoscillatory responses . 
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a l ' 

P 1 

a 

iwr 2 

2.20 

.25 

.969 

.1135 

1.69 

.3 

.931 

.2022 

1.238 

.35 

.894 

.236 

.811 

.4 

.859 

.2147 

.401 

.45 

.826 

.13625 


TABLE 4.1 


Example of Range of w(0) j | for Exponential Behavior 
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4. 8 Summary 

This chapter has presented a collection of analyses: of the J5MAC 
method. It has relied heavily on the special cases given in Section 3,1 
and as such, any conclusions with, respect to more general cases must Be. 
viewed with some caution. However, these special cases have Been care^ 
fully selected to accentuate types of behavior observed in more general 
settings and so it is believed that the qualitative results we have 
obtained are of some general interest. Thus, the result of our analysis- 
has yielded insight into the nature and causes of the behavior of the 
MMAC method. This understanding can be used to provide a basis for 
improving the design and for determining when it will work well. The 
main conclusions of this chapter are: 

1) At best the MMAC system is neutrally stable about an equilibrium 
point in that the probability has no tendency to return to its initial 
value following a perturbation (see Section 4.1) , 

2) If A CP) is an unstable matrix for P = l/2, then for the structure 
of Section 3.1 small perturbations result in the probability oscillating 
which in turn results in behavior which is either stable or unstable. 

If A CP) is a stable matrix for P = 1/2 either oscillatory or nonoscillatory 
behavior may occur depending on the size of perturbations and the 
stability of A(p) for P=0 and 1 Csee Section 4.1) . A procedure has 
been presented in Section 4,7 for determining a lower bound on the set 
of perturbations which yield oscillation-free responses, This: procedure 
is valid for any two model structure with- slight modification. 
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2 

3) The rate of change of P is proportional to | [w | | while the 
rate of change of | jwj | is proportional to | [wj j to the first power 
{see Section 4.1). This results in P changing faster than | |sj | for 

| |w| | large and slower for | |w| | small. This partly causes the square- 
wave behavior often noted for P. 

4) A necessary condition for | jwj | to converge to zero 

shown to be that A(P) be a stable matrix for some value of P (see Section 
4.1, 4.6) . 

5) A sufficient condition for exponential convergence has been 

/ 

shown to be that A(P) be a stable matrix for all values of P (see Section 
4.3). This, of course, is very restrictive. Lacking this condition, the 
possibility of at least short-term oscillations must be recognized. Thus 
the MMAC method is probably a poor choice when such oscillations can not 
be tolerated. 

6) If one model matches the true system and each model is diagonal, 

then | [w_j | 0 exponential (see Section 4.3).. This says nothing about 

the behavior of P(*) however. 

7) Whenever the MMAC method is used with two models, the equation 
for the probability can be divided into a scalar, static nonlinearity 
and a summation (i.e., the log likelihood ratio (see Section 4.4). Al- 
though not done here , this can be generalized m the N-model case to 

an N-l variable static nonlinearity and N-l log likelihood ratios. This 
approach is important as it emphasizes the switching behavior of the proba- 
bility, allows relatively siiple analysis to be done for the often seen 
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case when P is nearly piece-wise constant and provides for a slight 
numerical superiority in applications (see Section 5.3}. 

8) For two special cases specific conditions for the stability or 
instability of the oscillatory mode have been presented (see Sections 
4.5 and 4.6) . These results agree well with simulations and with each 
other. The basic result is that stability results when the stable and 
unstable eigenvalues are such that ab < 1 where a and b are appropriately 
defined as m Sections 4.5 and 4.6. The qualitative conclusion which 
follows is that for stability to occur the most unstable mode must be 
dominated by the stable modes. 

9) The oscillatory response may be stable even if no value of P 
results in A(P) being a stable matrix (see Section 4.6) . In this 
case the phase-plane plots resemble hyperbolas (see Figure 3.5) . This 
emphasizes the nature of the controller in that it attempts to achieve 
stability by alternately controlling each mode of the true system. 

10) For the special case of Section 3.1 with g = 0, it has been shown 
that v ( - ) is decreasing for P=0 and 1 as long as the LQ design procedure is 
used with a nonzero state weighting penalty. Generalizations have 
evaded analysis but this is believed to be a result of the gain and 
phase margin properties of the LO design [5]. 

Recall that in this special case the KF dynamics do not enter 
into the closed-loop behavior (since g = 0} . Thus, since the gam 
and phase margin properties can disappear when a KF is included in 
the feedback path [33], the extension to the general case (g^O) is 


in doubt. 
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11) In contrast to the previous conclusion, there are many values 
of the control gain which yield satisfactory matched behavior (i.e., 
stable') but which will result in unstable behavior in the MMAC system. 
Thus , care must be exercised in the choice of control gain. 

Thus, two major types of behavior have been observed and analyzed; 
oscillatory and exponential (nonoscillatory) . It has been seen that 
oscillatory behavior is very natural for the MMAC algorithm and that the 
conditions for excluding it are restrictive. Stability conditions for 
special cases have been presented and the qualitative impli cations for 
more general systems discussed. Thus, although based on special cases, 
the qualitative conclusions regarding the types of behavior and their 
underlying causes are believed to be of fairly general applicability. 



CHAPTER 5 


NUMERICAL ASPECTS OF MMAC 


In the course of this study, various issues relating to the numerical 
sensitivity of the MMAC method when implemented on a digital computer have 
come to light. The purpose of the present chapter is to discuss these 
aspects and finally to propose an alternate formulation which appears to 
overcome some of the limitations. 

In Chapter 4, while discussing the oscillatory behavior so often 
observed, it was noted that a lower limit on the probability tends to 
modify the behavior to some degree. Thus, while the analysis of Section 
4.5 indicates that the peaks of the state trajectories should be decreasing 
for the simulation in Figure 3.2, in fact, the peaks are seen to be 
constant. 

In Section 5.1 the oscillatory behavior is again analyzed, this time 
assuming that a lower limit on p is in use. This leads to the conclusion 
that if the limit is consistently achieved, then the peaks of the state 
trajectories will be equal when they would otherwise be decreasing. 
Furthermore , complete characterizations of 'both the period and peak 
amplitude of the states are given for both stable and unstable operation. 
The results are, except for the peak value, essentially the same as 
those in Section 4.5. 

Section 5.2 contains a brief discussion on various forms of the 
probability equation. It is found that a few forms behave very poorly 
in the face of numerical roundoff caused by finite precision. Precision 
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ref ers to the number of digits in the mantissa of the representation of 
a number in the computer. In addition, the need for extended range 
is also discussed. Range refers to the available size of the exponents 
in the computer. In particular, it is noted that due to large changes 
in the probability during oscillatory operation, (as seen in the simulations 
of Case 1 in Chapter 3) exponents of less than -75 often occur. Through- 
out this section the comments are made assuming that the algorithm is 
implemented using floating-point number representations. 

In light of the results contained m Sections 5.1 and 5.2, it is 
clear that the MMAC method places large demands on the computer which 
calculates the control. However, as seen in Section 4.4, there is an 
alternate form to the probability equation. In Section 5.3, a discussion 
of the numerical properties of this form is given which concludes that 
using Equation (4.11) in order to calculate the probability results in 
reduced sensitivity to the problems discussed in Sections 5.1 and 5.2 
compared to using Equation (2.22b). It is important to note that the 
proposed formulation does not eliminate the necessity of having a 
lower bound on the probability but merely allows the designer much 
greater latitude in the choice of this parameter. 

One final comment is important. The remarks in Section 5.2 are 
not based on any detailed analysis of the numerical 'properties . The 
purpose of the thesis has been to examine the stability of the MMAC 
algorithm and not to prove anything about its numerical properties. 

However, m performing, the simulations contained m Section 3.2, 
certain aspects of the numerical problem have become apparent. 
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The purpose of Section 5.2 is to point out these aspects. In keeping 
with this , no comment is made as to the numerical properties of the Kalman 
Filter, which is a significant problem in its own right. See, for ex- 
ample, the work of Sripad [30] . 

5.1 The Lower Limit Effect 

The basic equations of the MMAC method have been presented in 
Chapter 2. As mentioned there, it is often necessary in practice to 
place an artificial lower limit on the probability of any model. 

This is usually accomplished by simply setting the probability 
equal to the lower limit whenever the probability is less than the 
limit. Inclusion of a FORTRAN statement such as 

IF (P (I) .LT. PLIM) P{I) = PLIM 

is an example of the use of such a limit. The principal reason for 
such a limit is that without such a device, the probability of any model 
could be come exactly zero due to round-off errors. After this occurs, 
the probability remains zero as shown m Section 4.1. This is detri- 
mental in cases in which the true model is changing and also in cases 
in which an oscillatory behavior occurs (such as in Section 3.2.2), 
for then roundoff may occur before the probability switches. In fact, 
the Case 1 simulations of Section 3.2 are ]ust such a case. It should 
be recognized that although the lower limit is somewhat artificial, it 
is required in some form whenever a computer with finite precision is 


used for control calculations. 
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As mentioned in Se.ction 4.5, the existance of a lower limit requires 
the modification of the analysis of the behavior of the method. The 
remainder of' this section contains the modified analysis . As in Section 

4.5, we study the canonical problem of Section 3.1. Further, we assume 
that the system is oscillating as in Figure 4.3 such that P is closely 
approximated by either 0 or 1. Note that this is consistent with Section 

4.5. 

Recall from Section 4.5 that the canonical problem (see Section 3.1) 
can be rewritten as 



The log likelihood ratio, a(*), can then be written as 

k 

Ci(k) = Z Z.2 ^^2^-2 (k ^ " 2-i 

i_l 


(4.8) 


(5.1) 


where $ is defined in Section 4.5. 

—i 

In Sections 4.4 and 4.5 it is argued that during oscillatory periods 
(i.e., when P(-) is alternately near zero and near one), then a(k) must 
behave as the solid line in Figure 5.1 and also that a(*) can be approxi- 
mated by a piecewise-exponential curve. It will now be argued that the 
affect of a lower limit on the probability is to change the trajectory 
of a(*) to that of the dotted line m Figure 5.1. 

In Section 4.4 a static relationship between P(') and a(*) has 


been derived as 



.25 




85008AW01 8 


r 'CL{ • ) without limit 


<2(- ) without limit 
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P(k) 


P(0) (3e' 

_ -j a(k) 

P(078e + (l-P(O)) 


(4.11) 


Thus, if the probability of each model is constrained to be always greater 
than some limit then it can be shown that this is entirely equivalent 

to constraining a(*) such that 


a ( • ) - a < a, . 

s In 


(5.2) 


where is the value of a( •) when P(-) = 1/2 (given by Equation (4.13)) 
and is given by 




(5.3) 


Thus, placing a lower limit on the probability can be seen to modify the 
plot of a( ♦ ) to the dotted line in Figure 5.1: Initially, until the limit 
is reached, the limited and non- limited behavior are identical, but when 
the limit is reached, further increases in the magnitude of a(*) are 
ignored until a( ♦) again decreases (in magnitude). Note that since the 
size of the term in Equation (5.1) which forces a(*) to return to zero 
is independent of the existence of a limit, a(*) must cross the axis 
earlier when a limit is in effect than otherwise. This of course will 


influence future periods. It is thus argued that during oscillatory 
behavior in which the lower limit on the probability is achieved, in 
analyzing the behavior of the MMAC algorithm the positive term of 
Equation (5 . 1) can be set equal to f 1 or P « 0 while the negative 

term can be set equal to -a, . for P ~ 1 (assuming p, . small) . 
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Thus, we see that for the analysis of Section 4.5, the effect of a 
lower Ixmit on the probability is to allow us to replace Equation (4.23) , 
repeated here for convenience. 



(4.23) 


by the equation 


a 


lira 



(5.4) 


where has been defined m Section 4.4 as the time of the probability 
transition from P = 1 to P = 0 , 


a. = max | J A. (P) | | 2 
1 P 

and (X. is the maximum eigenvalue of d>. (where, by symmetry, (J - 0 - C) . 

1 ■*! ±2 

Noting as in Chapter 4 that | |y^ (T^) j | ~ j |y 2 (0) j j , Equation (5.4) 

can be rewritten as 


aa„ 




(5.5) 


or, solving for | (T^) 


li2<V 


a. 


I 2 “ - lll 2 (0)|| 2 


(5.6) 


where 


* 2 ° 


h - (a 2 -» 


(5.7) 


In the common case (see the simulations of Case 1 in Section 3.2) in 


which 
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||y 2 (°)|| 2 « H^tVll 2 

then Equation (5.6) becomes 


(5.8) 


% (T i> I 


a, 


lim 


*1 


(5.9) 


For the Case 1 parameters, this indicates that the probability will change 
when J |y (T^) | j r 100 which agrees well with the simulations when the 
integer nature of is considered. Note that Equation (5.9) implies 
that the probability will change when | jy^ (0) j | equals a constant if 
the lower limit on the probability is achieved; this constant is dependent 
only on the system parameters. 

In order to determine any change in period due to changes in the 
state, consider the situation in Figure 5.2 m which P = 1 for k = 0 to 
T^ and P = 0 for k = to where and are defined as the 

switching times of the probability. Thus, from Equation (4.20) 


ajM|y (0)|l 2 k e (0,T ) 

I iyqCte) 1 1 2 £ k . T T 

a 2 1 ^ ilZi (0, M 2 k e [t^ + t.,) 

where , 

\ = ! 11* 1 1 2 < 1 

a 2 = I|a 1 C0)|| 2 > 1 

by the construction of A (P) . From Equation (5.9) it is clear that 

!|y. 1 <o )|| 2 = IIz.i(t 1 + t 2 )II 2 or 
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T T. 


| 1^(0) | | 2 _< 3^3^ i 1^(0) 


(5.10) 


which-,- after -some manipulation , yields 



In 

In a 2 


(5.11) 


Evaluation of this results in the following conclusion: If a^a^ < then 
the peaks of the trajectories of | ]y^(«) | | are constant and the period 
1S increasing on each cycle of the probability. Furthermore , if = 1, 

then the peaks are constant and the period approaches a constant. The 
case in which > 1 is indeterminant from the above analysis although 

a plausible result, which is consistent with the simulations and the 
results of Section 4.6, is that the period would decrease until Equation 
(5.8) would no longer be a reasonable approximation. At this point, 
the peak amplitude of [ |y (•) | | would grow. Basically m this case the 
states are destablized more than they are stablized and the overall 
system is unstable. Neither terra (i.e., the positive or negative term) 
m the equation for a(*) (5.1) need become small. Note that a^a^ 1 
results in stable behavior which is in agreement with the results m 
Chapters 3 and 4 . However , m this case , note that the peaks remain 
constant and so asymtotic stability does not result. In fact, the 
oscillations will continue forever with 'constant amplitude and increasing 
period. This is indeed an unusual type of behavior. 
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5.2 Forms of the Equations and Accuracy 

In the course of doing the simulations contained in this thesis, 
a few points regarding the numerical aspects of the MMAC method have 
become apparent. These effects will not be dwelt upon as whole volumes 
could be written. However, they do form part of the basis for the form 
of implementation advocated in Section 5.3 and since such issues are 
often overlooked, a brief discussion is in order. 

For the simulations in Chapter 3, the variable P makes large 
excursions, in some cases changing by 5(3 orders of magnitude in one 
or two time steps. This, of course, makes implementation very crucial. 
For example, throughout this thesis, the identity 


P 2 = 1-P 1 (5.12) 

has been used to reduce the number of variables under consideration. 
However, using Equation (5.12) in a recursive equation such as 


P 1 (kjptr^) 

V k+1) = P 1 Dc)p(r 1 ) + (l-P^kHp'Cr^ 


(2.22b) 


can result in very poor accuracy due to roundoff. Thus, on a computer 
with 6-digit accuracy, this roundoff will occur for P t* ) ~ 10 re- 
sulting in P 2 ( * ) becoming exactly zero and remaining there unless some 
form of lower limit is placed on the probabilities as discussed in 
Section 5.1. Thus, using P 2 (0 = (1-P^*)) places severe restrictions 
on the choice of a lower limit on the probability, which may be un- 
desirable. The following three conclusions then are apparent. 
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1) Use of separate equations for P (.•) and P (•) is desirable if 

i. 

a recursion on P is to be performed as in Equation (2.22b). Thus, m ( 

pla ce of Equation (2_._22b) it is preferred, to use the pair of equations 

P^ (k)p (r^) 

P 1 (k+1) - + P 2 (k)p(r 2 ) 

P 2 (k)p(r 2 ) 

P 2 (k+i) - + P 2 <k)p(r 2 ) 

which avoids the direct recursive calculation of [l-P^(k)]. Thus, 
when using a floating point representation with finite precision ? 2 (k) 
in the above may be nonzero even when [1 - (k) ] is zero. Note that this 

does not significantly increase the computational load since both the 
numerator and denominator of the second equation are contained in the 
first. An alternate approach to solving this problem is discussed in 
the next section. 

2) Due to the fact that large variations' in both the probability and 
the states estimates often occur, the maximum number of digits possible should 
be maintained. It should be noted that this is very important for the 
probabilities but less so for the other variables since the rounding 

off of the state does not adversely affect the overall system. 

3) There is a direct relationship between the smallest possible 
value for p ^ and the possible range of exponents in the floating 
point representation. This tradeoff between design freedom (i.e., 
choice of P_ ) and computational complexity (i.e., needing large 

LXrl 
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magnitude exponents) must be recognized m any application. 

As shown in the first section of this chapter, the achieving of 
a lower limit by the probability can result in some change m the 
behavior observed in a system. Also, a number of the precision aspects 
of the problem have been shown to most directly affect the probability 
variable. It thus is of some advantage to develop an implementation 
which minimizes the effects of such phenomenon. Such an implementation 
is discussed in the next section. 

5.3 Proposed Implementation 

As shown in the previous sections of this chapter, the form of 
the equations used to implement the MMAC algorithm can affect the achieved 
results. The major area in which problems occur has been seen to be 
in the calculation of the probability. In this section an alternate 
form of the probability equation is proposed which, although not elimi- 
nating the problems, does tend to minimize their effect. 

The existence of a lower limit on the probability, although somewhat 
artificially imposed, is a necessary result of using a computer with 
finite precision. Thus, in any application of the MMAC method, a lower 
limit is present, in one form or another. The form of the probability 
equation about to be presented does not eliminate the lower bound problem 
but instead reduces the smallest possible value beyond that possible 
with the form given in Chapter 2 . It thus allows increased freedom of 
choice for the designer. 

From Section 5.2, one possible implementation of the probability 


calculation is 
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P.(k) P(r ) 

P±(k+1) p^dcjp^) + P 2 (k)p(r 2 ) 


i = 1,2 


It is argued there that this form is to be preferred to the form given 


by Equation (2.22b) in that the subtraction 1-P^(*) is avoided. In 
Section 4.4, it is shown that Equation (2.22b) is entirely equivalent 
to 


P x (k) = 


P 1 (0) Be 


- fa (k) 


-|a(k) 

P 1 (0)Be +P 2 (0) 


(4.11) 


with 

k -1 -1 

a(k) = 2 lr^(i)9 1 i r l (i) - r' (ije^rd)] . (5.13) 

i=l 

It is believed that Equation (4.11) is numerically superior to either 
of the other approaches. Note that, while a term such as Il-P^(O)] is 
still required to compute ? 2 ( - ) for the calculation of the control, any 
errors in this computation do not get accumulated in the probability - 
Thus, the real advantage of using Equation (4.11) is that it is static, 
involving no recursive calculations. Thus, while roundoff obviously 
will occur, it will not influence future values as would happen with the 
other approaches. It should be pointed out that the recursive nature 
of the probability has, in effect been retained in the calculation of 
a(*) and that now care must be taken to guarantee its numerical accuracy. 

However, by its nature, it tends to be a better-behaved function. For 

-50 

example, a change in P from 1 to 10 corresponds to a change m a(*) 
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from approximately +225 to -225. Thus, the recursion has effectively 
been moved to the exponent (i.e., the Log Likelihood Ratio) rather than 
in the probability itself. The importance of this can best be seen by 
considering the following problem, assume one wishes to calculate the 
quantity (e^e^) . Once approach is to compute (e a ) , then compute (e* 5 ) 

cl J) 

and finally (e ) (e ) . This is effectively the kind of approach Equation 
(4.22b) uses. However, it is more accurate (for finite precision calcu- 
lations) to compute e^ a+ ^^, which is analogous to what Equation (4.11) 
does. 

•A further advantage of using Equation (4.11) is that it can be approxi- 
mated by any one of a number of functions without having to be overly con- 
cerned about the accuracy of the approximation. This is due to the fact 
that Equation (4-11) itself is not recursive and thus errors made in ap- 
proximating P do not accumulate except through the true states. One ex- 
ample of such an approximation is a switch such that P ( • ) = 0 for a ( • ) >0 and 
P = 1 for Ct ( * ) < 0 while another has been presented in Section 4.5. 

Note that for the variations of the MMAC algorithm which use the control 
for the most probable model rather than the weighted control, a switch 
is an exact representation and not an approximation. 

It should be pointed out that this approach can be generalized to an 
N model situation. In this case, one can define the N-l log likelihood 
ratios based bn, for example, model N as 

k -1 -1 
a. (k) = S [r'(j)0 r C j ) - r' tj ) 9 £ N ( 3 ) J 

1 . , — i — 1 —1 — N — N 

3=1 
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where all terms are defined analogously to their two -model- case counter- 
parts. Repeating the derivation in Section 4.4 yields the N-l nonlinear 
equations for the probabilities 


P ± (k) 


P. (0)3. e 

i i 


K (k > 


N-I 

s 

3~1 


P (0)3. e 
3 3 


- (k)l 


+ El-P„(0)] 

N 


5.4 Summary 

In this chapter a number of issues relating to the behavior of the 
I4MAC algorithm when implemented on a digital computer have been examined. 
The most important of these is contained in Section 5.1 m which the 
effect of placing a lower bound on the probability is discussed. The 
principle conclusion of that section is that during oscillatory periods 
in which the lower limit is achieved, the conclusions of Chapter 4 are 
still valid except that the peaks of the trajectories of | jy^j | are con- 
stant even during stable operation. Table 5.1 thus summarizes the con- 
clusions about the stability whenever the lower limit is achieved. These 
results are seen to agree well with the simulation studies in Section 
3.2.2. 

It should be born in mind that, as discussed earlier in the chapter, 
there are philosophical reasons for choosing a relatively large value 
for the variable P„„. Note that such a choice results in a tendency to 
reduce the amplitude of oscillations as well as the period since the 
unstable state must now change less to cause a ( • ) to change by a T/7? j 
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(i-e., cause P to reach P T J. 

LIM 

Section. 5.2 has presented a brief discussion of the numerical 
properties of various forms of the probability equation. The principle 
conclusion is that the MMAC algorithm, if implemented using the straight 
forward equations such as Equation (2.22b), places large demands on a 
digital computer due to the recursive nature of the calculations , the 
use of a term such as (l-P(k)) at each iteration, and a possibly large 
dynamic range of the probability. An alternate formulation, based on 
the analysis of Section 4.4 is presented in Section 5.3. It is shown 
that this formulation tends to be superior to the others presented in 
that the recursion is done on the exponent where accuracy can be better 
controlled. 


v 



CHAPTER 6 


COMMENTS AND CONCLUSIONS 

Tins thesis has presented an analysis of the behavior of the Multiple 
Model Adaptive Control algorithm through both analysis and simulation. 
Extensive use has been made of a canonical system (see Section 3.1) in 
which certain assumptions on the structure of the true model are made, 
chief of which is that it is diagonal. While somewhat extreme, this 
sample structure has been carefully chosen to display what we feel are 
the critical characteristics of the method as observed m more general 
applications and which, unlike the more general problem, also has been 
amenable to detailed analysis. 

The purpose of this chapter is twofold. First of all, various 
ad hoc modifications have been proposed for the MMAC algorithm which 
are aimed at overcoming some of the undesirable behavior observed in 
applications. An example of such a modification is the introduction of 
a low pass filter to smooth out the probabilities. Section 6-1 contains 
a brief discussion of the most prominent modifications along with three 
methods for improving the response which are suggested by the results 
of the analysis of this thesis. The basic conclusion is that with any 
of these modifications great care needs to be taken to ensure that the 
response is not degraded. 

The second purpose of this chapter is to provide a summary of the 
major conclusions of this thesis. This is done in Section 6-2 m which 
the specific conclusions of the analysis contained in this thesis are 
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summarized. Also in this section are a list of qualitative conclusions 
regarding the MMAC system when applied to general systems. These result 
from extrapolating the specific conclusions to more general systems. 

The chapter concludes with a list of suggestions for future research in 
understanding the MMAC algorithm. 

6.1 Modifications 

Various ad hoc modifications have been proposed -to overcome some of 
the undesirable properties of the MMAC algorithm which have been observed 
in applications. In this section two of the most prominent of these are 
discussed. First of all, in order to make the algorithm more sensitive 
to changes in the true model, a form of the MMAC algorithm which possesses 
a finite memory property can be used. In the first subsection, it is 
shown that most of the properties of the MMAC method apply in this case. 
The exception to this is that no_ convergence property can be given due 
to the finite memory. 

The next subsection contains a brief discussion of the addition of 
a low pass filter in the probability calculation. This modification 
has been proposed [23j to smooth out the rapid probability transitions 
and oscillations which can occur and thereby attempt to get a smoother 
state response. It is argued in Section 6.1.2 that, especially when no 
model matches the true system, this can result in at best no change in 
performance and at worst in a destablilizing effect. This section then 
concludes with a discussion of three modifications which have been 
suggested by the analysis of this thesis. 
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6.1.1 Finite Memory MMAC 

One of the many ad hoc modifications to the MMAC method which has 
been proposed has been termed Finite Memory MMAC (FM-MMAC) . It can best 
be understood by considering the form of the probability equation given 
m Section 4.4: 


P(k) 


« - 7 ct(k) 

P(0)6e 

-- a(k) 

P(0)ge + (1-P(0)) 


(4.11) 


where 

a(k) = S r^(i)0. 1 ' L r 1 (i) - (i) (i) • (4.10) 

i_l 

This form of the equation makes it clear that all of the past data are 
equally weighted in determining the present probability. This is not 
necessarily what one would desire in an adaptive controller since data 
from the far past may not be relevant to the current operating conditions. 
There are, of course, many ways to modify the method so as to make it 
more responsive to the immediately past data than" it is to the more 
remote data. For example, one may exponentially weight the past data: 

a(k) = S a (k_l) [r^(i)0^E (i) - r’dJ^d)] (6.1) 

i_l 

for some a < 1. An alternative is to merely use the last M values of 
the residuals. Thus, one could use 


5(k) = 2 [r^ije^^Ci) - r’ (i)6 2 1 r 2 (i)3 

i= k-M+1 


( 6 . 2 ) 
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where M is a free design parameter which would be chosen partly on the 
basis of how fast the parameters of the true model are expected to 
change. This is referred to as the Finite Memory MMAC algorithm. 

It should be noted that a significant amount of the analysis of 
systems of these types has already been done m Chapter 4. To illustrate 
this, a discussion of the second of the proposed schemes with M=1 is 
given. Thus, take 

a(k) = r^ (k) (k) “ £.2 2 - 2 ^^ * (6.3) 

First of all, from Section 4.1, the basic linearized analysis con- 
clusions such as the neutral stability are unchanged as are the results 
of Sections 4.2 (Universal Stability). Note that the state convergence 
results of Section 4.3 do not carry over in this case. The technical 
reason is that ( • ) is no longer non- decreasing. In fact, it is easy 
to show that 

i |w(-) [ j 0 => P ^ 

where N is the number of models. Thus, if the closed loop system is not 
stable for the case in which all models are equally likely (i.e., the 
N-model extension of the condition that A(P) is stable for P = 1/2) , 
then | | w | | can not approach zero . This is a direct result of the finite 
memory assumption and occurs for any finite M. This is not , however, 
necessarily bad for an adaptive controller. It roughly corresponds to 
the notion that without information as to which model is correct. 


(i.e. , 


r. = 0) one can assume that each model is equally likely. 
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The results of Section 4.4 on the analysis of the probability 
equation also hold since they do not depend on how a ( * ) is computed. 
Thus, Figure 4.1 is again useful as are the approximations 


a(') 

» 

a 

s 

=> 

P(*) 

s 0 

<*(•) 

« 

a 

s 

=> 

P(-) 

~ 1 

and (4.14), where 

a 

s 

is as 

in Equation 


(4.12) 


can best be appreciated by noting that combining Approximation (4.12) with 
Equation (6.3) yields 


Ir^k)!! _ x » ||r 2 (k)|l =>P(k) * 0 

% ®2 


1^00 I I « ||£ 2 (k)|i =>P(k) * 1 

£l ®2 


(6.4) 


where, as in Chapter 4 we have assumed => 0. (Removing the assumption 

that a - 0 merely results m a being added to the appropriate term in 
s s 

Approximation- (6.4) . ) Thus one can expect that, compared to the infinite 
memory case, the period of oscillation for the M=1 case will be shorter 
due to a lack of the accumulation of the residuals over time. Effectively, 
by removing the summation of the residuals (the discrete time analogue 
of an integrator) , a lag has been removed from the system. Thus, one 
would expect the peaks of the state curves to be smaller than when the 
summation is included. Note that it is possible, following the analysis 
of chapter 4, to calculate approximations to the switching times T^ and 
T 2 - In fact, it becomes somewhat easier to compute since no summation 


is involved, Note, however, that care must be used in applying the results 
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due to the assumption that P is always near zero or one. 

Simulation results of this case are shown in Figure 6.1. The con- 
ditions of this simulation . are identical to those of Figure. 3 . 2 .except. . .. 
that M=1 is used in the probability calculation. Note, however, that 
the character of the response has changed markedly. The probability 
(Figure 6.1a) changes much more rapidly and takes on intermediate values 
occasionally as it transitions between zero and one. This, of course, 
complicates any analysis such as the computation of the switching times 
and . The state trajectories are shown in Figure 6.1b. Note that 
they are also much more oscillatory and in fact appear to limit-cycle. 

It is unknown at present if this is in fact a true limit cycle. Also, 
clearly shown is the reduction in the peaks of the state trajectories 
from about 17 for the unmodified situation to about 1.2 for the M=1 
case. 

The variable In (x x„ ) , shown in Figure 6.1c, also is highly 

-i* ^ 

oscillatory and does not exhibit rhe negative slope which is evident 
in the unmodified case (Figure 3.2). This is further evidence of a 
limit cycle and indicates that the closed loop system is only neutrally 
hyperbolically stable. It should be noted that the analysis of’ Sections 
4.5 and 4.6 on the conditions for hyperbolic stability required that 
the probability be always near zero or one. As long as this condition 
is met, the results of those sections still apply to the Finite Memory 
case. Note, however, that the assumption on the probability may not 
be as reasonable due to the presence of the high frequency oscillations. 
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6.1.2 Low Pass Filter 

As seen m Section 3.2, the probability often exhibits a zero-one 
type behavior, rapidly changing between these values. In an attempt to 
smooth out these probability transitions, it has been proposed that a low 
pass filtered version of the probability be used for the control calcu- 
lation. In this section a brief argument is presented which concludes 
that this may result in a degrading of the response. 

As seen in Chapter 4, when A{P) is an unstable matrix for all fixed 
P, the oscillatory behavior is necessary if hyperbolic stability is 
to be achieved. Thus, one can easily imagine cases in which low pass 
filtering the probability would result in a degradation in response. 
Although no detailed analysis has been done, the results of this thesis 
indicate that at least during the initial period, including a low pass 
filter on the probability would cause the probability to, lag behind 
what it would be otherwise. This in turn would result in the unstable 
state growing to a larger value before the probability switches. Although 
the effects in later periods are unknown, it is reasonable to conclude 
that the net result is an increase in the values of the peaks of the 
state trajectories compared to the case in which no extra filtering is 
done. 

6.1.3 Miscellaneous Modifications 

In this section three relatively simple modifications to the MMAC 
algorithm which have been suggested by the analysis contained in Chapters 
4 and 5 are introduced and discussed. The three modifications are: 
scale the residuals, increase the lower limit on the probability and take 
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the square root of a ( ♦ ) in the probability equation. 

In Chapter 4 it is shown that a discrepancy between the probability 
equation and the state equation exists in that for P near 1/2 the rate 
of change of the probability is proportional to the square of [ jwj | while 
the rate of change of | |w| j is only proportional to | jw| j to the first 
power. This suggests that one way of smoothing the probabilities is 
to replace Equation (2.22b) by 


P L (k+1) 


P 1 (k)p(r i ) 

P 1 (k)p(r i ) + ? 2 (k)p (r^) 


(6.5) 


where 

- ~ Vr : eT^r 

p(r ± ) = ^e 1 r 1 . (6.6) 

Note that p(*) is no longer the Gaussian density function. This version 
will be referred to as the square root modification of the MMAC algorithm. 

Simulations using this redefined function (shown m Figure 6.2) 
indicate that the probability still behaves in a zero-one fashion. 
Furthermore, note that the peaks of the state trajectories have increased 
significantly compared to the results using the normal probability 
function. (See Figure 3.2). This is believed to be due to the effective 
delay which results from C4( • ) changing linearly rather than quadraticly. 
Thus, this modification does not accomplish the goal of smoothing the 
probability and results in degraded state responses. It is further 
concluded that the dominant nonlinearity m so far as the oscillatory 
behavior is concerned is that of Equation (4.11) 
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P(k) 


- \ a(k) 

P (0) Be 

- \ a(k) 

P(0)Be + (1-P Q (0)) 


(4.11) 


and not those of the a(*) term. 

Finally, note that as predicted by the results of Chapter 4, the 
closed loop system is hyperbolically stable. This is best seen by con- 
sidering the plot of InCx^x^) given in Figure 6.2c in which this variable 
is seen to decrease linearly. This indicates that the states behave in 
such a fashion that 


* 1*2 = 6 


bt 


where b is the slope of the line in Figure 6.2c. 

The second possible modification is to change the lower limit on 

the probability, P, , discussed in Section 5.1. In that section ir. 

lim 

is shown that during periods when the lower limit is achieved, then the 
peaks of the state trajectories are given by where 


d, - 2 In 
Inn 


1-P. . 
lxm 

P, 

lim 


(5.3) 


and is given by Equation (5.7). Thus it is possible to reduce the 
peaks merely by increasing the value of P ^ a s long as the approxima- 
tions 


A (1) « A (1-P_ . ) 
— — lim 

A (0) « A (P ) 
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are valid. Even when these approximations are not valid, the approach 

is still valid if the effect on R^ is accounted for by replacing 

I IS. (!) | | with | jS, (1-P, , ) | | (see Section 5.1) . 

— i ' ' 1 '—i lna 1 1 

Simulations m which is changed have been done and they confirm 

this observation. However, two points need to be recognized in exploiting 
this fact. First of all, if P re ^ uce< ^ too much, the assumption 

that a ^/R is greater than | |y 2 (0) j| 2 (see Section 5.1) is violated 
resulting in the peaks not being as predicted by a m / R 1 . Also, by 
Equation (5.3), the value of the peak varies as the natural logarithm of 
P lim SO large changes m P are necessary to make moderate changes 

in the peaks. 

A third modification, suggested by the analysis of Section 4.4, 
is to scale the residuals . As seen in Section 4.4, the value of the 
probability is determined by a ( • ) . Thus, if the calculation of a(*) is 
replaced by 

&(*) = a a(*) (6.7) 

with a>l where a(*) is given by Equation (4.10), then the changes in 
at*) are amplified. This, when combined with the existance of a lower 
limit tends to reduce the peak value of the state trajectoreis and 
the period (see Section 5.1). A simulation of this, shown in Figure 
6.3, indicates that significant decreases in the peak amplitudes are 
possible. The simulation conditions are identical to those of Figure 3.2 
except that the residuals are scaled by a factor of 10. This corresponds 
to letting a = 100. in Equation (6.7). As a comparison, the unmodified 
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MMAC algorithm had a peak of about 17. (see Figure 3- 2b) , the FM-MM&C 
algorithm had a peak of about 1.4 (see Figure 6.1b) but oscillated rapidly , 
while scaling with a value a = 100 results in peaks of about 1.6 (see 
Figure 6.3b). Note also that the quantity In x^x^ for the case when 
a(* ) is scaled (Figure 6.3c) behaves exactly as in the unmodified case 
while the FM-MMAC algorithm, results m neutral hyperbolic stability 
(i.e. , m In x^x^ being constant) . 

This modification only reduces the peaks of the state trajectories 
when the lower limit on the probability is achieved. In fact, it follows 
directly from the anslysis of Section 5.1 and thus all of the analysis 
of that section can be repeated. In particular, the results on the 
hyperbolic stability and the calculation of the switching times follows 
analogously with the use of Equation (6.7). One drawback of this 
approach is that since a(*) is scaled without modifying the associated 
KF's, a(* ) will tend to be affected more by noisy inputs than before. 

This in turn could significantly alter the response to noise m the 
observations (see Section 6.3). 

6. 2 Summary of Results 

The major results of this thesis are of two basic types. The 
first type is the specific conclusions which have been derived for the 
special cases discussed in Chapter 3 and 4. These result in specific 
criteria for each type of behavior but, strictly speaking, apply only 
to the associated special case. Of* possibly greater importance are 
the qualitative results which, while not proved in this thesis are 
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reasonable extrapolations of the specific conclusions. Both types of 
results are summarized below. For the details of the specific conclusions, 
the reader is referred to the appropriate section of the body of the thesis. 

The major specific conclusions of this thesis are: 

1. At best the MMAC system is neutrally stable about an equilibrium 
point in that the probability has no tendency to return to its initial 
value following a perturbation. (See Section 4.1). 

2. If A(P) is an unstable matrix for P=l/2, then for the structure 
of Section 3.1 small perturbations result in the probability oscillating 
which m turn results in behavior which is either bounded or unstable. 

If A (P) is a stable matrix for P = 1/2 either oscillatory or nonos dilatory 
behavior may occur depending on the size of the perturbations and the 
stability of A(P) for P- 0 and 1 (see Section 4.1). A procedure has 
been presented in Section 4.7 for determining a lower bound on the set 
of perturbations which yield oscillation- free responses. This procedure 
is valid for any two model structure with slight modification. 

3. The rate of change of P is proportional to j |w[ | while the 
rate of change of | jw| j is proportional' to | |w;( j to the first power (see 
Section 4.1). This results in P changing faster than | |w| | for j |wj | 
large and slower for- | [ w | | small. This partly causes the square-wave 
behavior often noted for p. 

4. A necessary condition for j |w| | 0 has been shown to be that 

A(P) be a stable matrix for some value of P (see Section 4.1, 4.6). 

5. A sufficient condition for the exponential convergence of the 
state has been shown to be that A(P) be a stable matrix for all values 
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of p. (see Section 4.3). This, of course is very restrictive. Lacking 
this condition, the possibility of at least short term oscillations 
in the probability which cause the states to increase must be recognized. 
Thus, the MMAC method is probably a poor choice when such oscillations 
can not be tolerated. 

6'. If one model matches the true system and each model is diagonal, 
then | |w| | -> 0 (see Section 4.3) . This says nothing about the behavior 
of P ( • ) however. 

7. Whenever the MMAC method is used with two models, the equation 
for the probability can be divided into a scalar, static nonlinearity and 
a summation (i.e., the log likelihood ratio) (see Section 4.4). Although 
not done here, this can be generalized in the N-model case to an N-l 
variable static nonlinearity and N-l log likelihood ratios. This ap- 
proach is important as it emphasizes the switching behavior of the prob- 
ability, allows relatively simple analysis to be done for the often 

seen case when P is nearly piece-wise constant and provides a slight 
numerical superiority m applications (see Section 5.3). 

8. For two special cases specific conditions for the stability 

or instability of the oscillatory mode have been presented (see Sections 
4.5 and 4.6). These results agree well with simulations and with 
each other. The basic result is that stability results when the slowest 
decaying stable eigenvalue and the most rapidly growing eigenvalue 
are such that their product- is less than unity. 

9. The oscillatory mode of response may be hyperbolicly stable 

m that the quantity ln(x^x^) is decreasing even if A(P) is an unstable 
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matrix for all constant values of P (see Section 4.6). This is seen to 
correspond to the observation that the states are alternately controlled 
in a manner such that 

bt 

X 1 X 2 = e 

where b is a constant which is negative for hyperbolicly stable systems. 

In this case the phase-plane plots resemble hyperbolas (see Figure 3.5). 
This illustrates the nature of the controller in that it attempts to 
achieve stability by alternately controlling each mode of the true system. 

10. For the special case of Section 3.1 with g = 0, it has been 
shown that v ( ■ ) is decreasing for P=0 and 1 as long as the LQ design 
procedure is used with a nonzero state weighting penalty. 

11. In contrast to the previous conclusion, there are many values 
of the control gain which yield satisfactory matched behavior but which 
will result m unstable behavior m the MMAC system. Thus, care must 
be exercised m the choice of control gain. 

12. The existence of a lower limit on the probability slightly 
affects the results. When the limit is consistently achieved, the 
peaks of the state trajectories tend to a precomputable constant (see 
Section 5.1). 

13. Various forms of the equation for the probability have been 
examined as to the effects of numerical roundoff (see Section 5.2) 
resulting m the proposal m Section 5.3 of a form which is believed to 
be superior to others in that it is less sensitive to roundoff effects. 
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14. Five ad hoc modifications to MMAC algorithms have been 
examined (see Section 6.1). The Finite Memory MMAC resulted in a 
significant reduction in the amplitude of the peaks of the state 
trajectories but also resulted in a significant increase m the fre- 
quency of the oscillation. Adding a low pass filter was seen to be 
possibily detrimental to the simulation results as was using the 
square root of a(*) in the probability calculation. Two approaches 
which appear to_ improve the response have been given. These -are to 
increase the smallest value that the probability is allowed to assume 

(i.e. # increase P, , ) and also to scale the residuals by a scalar which 
lun 

is greater than one. 

The qualitative results are: 

1. If one of the' models matches the true system, good behavior 
(as. defined by the design cost function) probably results. Some 
doubt must remain due to the fact that convex combinations of stabi- 
lizing controls do not necessarily yield stabilizing controls. 

2. If none of the models match the true system but if at least 
one model results in a stabilizing control (that is, if A (I?) is the 
extension of A(P) to the N model case, then this results when there 
exists some i such that P_^ = 1 => A (P) is a stable matrix) then 
stability probably results. However, it is clear that the quality of 
the response may still be poor in that it may be signifiantly different 
from the design performance specifications as measured by the response 
of the control system under perfectly matched conditions. Furthermore, 
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it is possible to get oscillations in this situation. Consider, for 

example, the case with two models where one model stabilizes the closed 

loop system but has a slow filter while the other destabilizes the system 

but has a fast filter. One could then expect a cycle to develop in which 

the system alternates between the stable and unstable control laws. 

\ 

3. Even if no constant value for the probability results in 

a stabilizing control, the overall system may still be hyperbolically 
stable (where hyperbolic stability would have to be redefined for the 
general case) . In this case, the probabilities would tend to rapidly 
change alternately controlling the various modes of .the true system. 

4. Even when a constant value of P results m a stable control, 
oscillations in the probability are likely to occur. In fact, the 
oscillations are a natural and necessary attribute of using the MMAC 
algorithm when perfect model matching does not occur. 

5. The MMAC algorithm can be expected to provide good response 
in conditions in which the set of possible models is closely approxi- 
mated by a finite set of models. An example of this would be the 
control of a system subject to discrete failures in various components. 

6.3 Suggestions for Future Research 

Many aspects of the behavior of the MMAC algorithm remain to be 
studied. A partial list follows. 

1. Ljung et al [24, 25 j have considered closed loop identification 
for some specific model structures and identification methods. Appli- 
cation of their results to the MMAC method remains to be done. This 
may result in the extension of Baram's results [15] to the closed loop 


case 
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2. In Section. 4.5, the analysis produced the conclusion that the 

peaks of the state trajectories decrease and the period increases if 

) 

a i a 2 1S ^ ess t ^ ian unity. However, this eventually results in the 
violation of the assumption that the probability is either zero or one. 
Furthermore, if the closed loop system with P constant is unstable for 
all values of P, then the state can not approach zero since, as seen in 
Chapter 4, this would imply that P approaches zero even faster than the 
state. Simulation results to date have been of little help in that in 
all simulations m which long term behavior was explored an underflow 
to zero occurred for at least one of the components of the state. Thus, 
the asymptotic properties of the MMAC method remain to be determined. 

3. All of the analysis and simulation results of this thesis are 
based on a noise-free case. The effect of noise needs to be determined. 
Note that this is a non-trivial task since even if all the r ( • ) are 

—i 

assumed zero mean, aC*) is then Chi-Squared resulting m a very complex 
distribution for the probability. 

4. Many of the qualitative properties discussed in Section 6.2, 
while substantially based on the specific conclusions of this thesis, 
remain to be proved in general. An example of this is the proof of the 
general convergence of the probabilities when one of the hypothesized 
models matches the true system but when the assumption that the' true 
system is diagonal is not satisfied. 

5. Much of the analysis contained m this thesis has been based on 
the canonical problem structure of Section 3.1 in which the true system 
is assumed to be diagonal with two states. Extensions to non-diagonal 
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systems and/or systems with more than two states are needed. 

6. Much of the analysis has included the assumption that there 
were two hypothesized models. Extensions to the N- mo del case are needed. 

7. In the analysis of Section 4.6 on hyperbolic stability it is 
assumed that the control gain corresponding to the mismatched state 

(g) is zero. Extensions to the more general case have been attempted 
without success (see Section 4.6). However, such results would yield 
useful insight into the general MMAC system. 

8. In Section 4.6 it is noted that for the special case under 
consideration, the LQ methodology is guaranteed to produce a control 
which meets one of the criteria for hyperbolic stability. It remains 
to determine if this is a general property or merely a result of the 
conditions of the special case. 

9. Recall from Section 4.6 that ~ hyperbolic stability is guaranteed 
when Conditions (4.41) and (4.42) are met. The first of these guarantees 
hyperbolic suability for P=0 and 1 while the latter does so for P = l/2. 

If it can be shown that in cases in which Condition (4.41) is met but 
(4.42) is not that P oscillates in such a fashion that it never takes 
on intermediate values, then Condition (4.41) would itself be sufficient 
to guarantee hyperbolic stability. 

10. In Section 6.1 it is seen that using the finite memory MMAC 
algorithm with M=1 results m what appears to be a limit cycle. Simula- 
tion results indicate that the characteristics of the limit cycle are 
dependent on the system parameters and not on the initial conditions. 
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For example, when initialized with extremely small initial conditions 
(i.e., true states of 1 x 10 i0 and 1 x 10 20 respectively) behavior 
very much like that in Figure 6.1 has been observed in that oscillations 
of the same amplitude and period start after about 200 time steps. It 
would thus be useful to examine the characteristics of this case to 
determine the amplitude, bias {see Figure 6.1) and frequency of this 
oscillation, 

11. Throughout this thesis, extensive use has been made of the 
matrix norm | }a| j = Max X(AA’). This results in conservative estimates, 
for example, for the domain of attraction in Section 4.7. Calculations 
in which Max|X(A) j is used m place of the norm have been done and 
estimates which agree better with simulation results thus obtained. How- 
ever, this quantity is not , in general, a norm for the matrix A. Research 
into using Max) A (A) | m place of the norm may result in the improvement 

in the estimates used throughout this thesis. 

12. In this thesis, the concept of hyperbolic stability is intro- 
duced for the two-state diagonal system. Extensions of this useful 

concept to more general system structures needs to be done. 

* 

13. A few design modifications have been proposed and analyzed. 

It is believed that the results of this thesis can be used to develop 
additional design modifications which lead to improvements of the 
response. 

14. The results of this thesis strongly indicate the types of 
behavior which a general MMAC algorithm can cause. Although no con- 
sistent general methodology has been developed, the qualitative results 



-167- 


of this thesis provide a basis for such methodology. Specifically, 
important design ground rules are needed for selecting an appropriate 
hypothesized model set given the expected range of parameters of the 
true model and for the selection of the control gains for each of the 
hypothesized models. Our analysis in Chapter 4 provides insights into 
the constraints which must be imposed upon such choices and thus should 
help form the basis for such a design methodology. 



APPENDIX A 


DEVELOPMENT OF THE ALGORITHM FOR THE RANGE 


OF 


W_ 


FOR EXPONENTIAL BEHAVIOR 


In Section 4.7 an algorithm has been presented which allows one to 
compute a bound on | Jw(0) j | such that exponential (non-os dilatory) 
behavior results . This appendix contains a more detailed development 
of that algorithm. 

Consider Equation (2.22) 


w(k) = A(P) w(k) 


P (k+1) 


P (k) p(r ) 

P (k) p (r^) + (1- p(k))p(r 2 ) 


We make the following assumptions as in Section 4.7: 


1) p(0} = 1/2 and 6=1 U.e, . a = Qi 

s 

2) A(P) is such that there exists an £ < 1/2 such, that 


1 |a(p) | ( < 1 VPS Cl/2 ^S, 1/2+ el 

j | A (P) | | >_ 1 otherwise . 

Define 6 to be the value of a ( * ) corresponding to P = 1/2- £; that is, 
define 6 as the solution of 


.5 - £ 


6e 






(A.l) 
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Thus, by the definition of e and a it can be shown Csee Section 4.4) 
that 


P(*) e (1/2 -e, 1/2 + £) <=>a(*) e (a + <5, a -a) (a. 2) 

s s 

Combining Equations (4.48) and (A. 2) with the definition of a(*) from 
Equation (4.10) , it can be shown that 


| | A(P (k) ) j | < 1 Vk > 0 

<=> 

00 

1 | ct (°°) | | = | 2 [r 1 (i)9^ 1 r 1 (x) 

i=l 


—2 (D— 2 1 — 2 (i^ ■* 1 < ® ’ 


(A. 3) 


This condition merely states that for A(P{»)) to remain a stable matrix 
P('), and therefore also a (• ) / must remain close to their initial values. 

It is now necessary to examine a ( • ) . Equation (4.10) can be rewritten 
as 


k 

a(k) = Ew' (0) 

i=l 


i 

n 

j=l 


A(P(i)) ’ 


l 


U A(P (i) ) 
D=i 


w(0) 


where 


0 


2 = 


0 


0 



Define a such that 


(A.4) 


|!a(P(i))|| 2 <a Vi = 1,2,... ,k 


Taking norms in Equation- (A.4) then yields 
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a(k) 


k 

Z 

i-1 



(A. 5) 


where a = MIN {A (9^) , A ( 0_^ ) } with A ( 6_) representing the set of all eigen- 
values of 0. Evluating the summation in Equation (A. 5) then yields 


a(k) 


< 



(1 - a) a 



(A. 6} 


Thus , the approach m the remainder of this section is to hound changes 
in a(*) so that A(P (*)) is always a stable matrix. If A(P{*))is always 
stable, then a is less than one and Inequality (A. 7) can be replaced by 

, , a liHoN 2 

|afr>! 1 (1 - a) cr Vk (A * 7) 

since the right hand side of Inequality (A. 7) is an increasing function 

2 

of k. This leads to the procedure of Section 4.7 for finding | |w^j j 
such that ja{*} j < 6 ¥k . 
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